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Abstract. It is well known that cosmic rays contribute significantly to the pressure of the interstellar medium in our 
own Galaxy, suggesting that they may play an important role in regulating star formation during the formation 
and evolution of galaxies. We here discuss a novel numerical treatment of the physics of cosmic rays and its 
^ ■ implementation in the parallel smoothed particle hydrodynamics code GADGET-2. In our methodology, the non- 

l/") ' thermal cosmic ray population of each gaseous fluid element is approximated by a simple power law spectrum 

in particle momentum, characterized by an amplitude, a cut-off, and a fixed slope. Adiabatic compression, and 
a number of physical source and sink terms are modelled which modify the cosmic ray pressure of each particle. 
The most important sources considered are injection by supernovae and diffusive shock acceleration, while the 
primary sinks are thermalization by Coulomb interactions, and catastrophic losses by hadronic interactions. We 
also include diffusion of cosmic rays. Using a number of test problems, we show that our scheme is numerically 
robust and efficient, allowing us to carry out the first cosmological structure formation simulations that self- 
• consistently account for cosmic ray physics. In simulations of isolated galaxies, we find that cosmic rays can 

significantly reduce the star formation efficiencies of small galaxies, with virial velocities below ~ 80kms -1 , an 
effect that becomes progressively stronger towards low mass scales. In cosmological simulations of the formation 
5_l ' of dwarf galaxies at high redshift, we find that the total mass-to-light ratio of small halos and the faint-end of 

^ | the luminosity function are strongly affected. The latter becomes flatter. When cosmic ray acceleration in shock 

waves is followed as well, we find that up to 40% of the energy dissipated at structure formation shocks can 
appear as cosmic ray pressure at redshifts around z ~ 3 — 6, but this fraction drops to ~ 10% at low redshifts 
when the shock distribution becomes increasingly dominated by lower Mach numbers. Despite this large cosmic 
ray energy content in the high-redshift intergalactic medium, the flux power spectrum of the Lyman-a forest is 
only affected at very small scales of k > 0.1km _1 s, and at a weak level of 5 — 15%. Within virialized objects, we 
find lower contributions of CR-pressure, due to the increased efficiency of loss processes at higher densities, the 
lower Mach numbers of shocks inside halos, and the softer adiabatic index of CRs, which disfavours them when a 
composite of thermal gas and cosmic rays is adiabatically compressed. The total energy in cosmic rays relative to 
the thermal energy within the virial radius drops from 20% for 10 12 /i _1 Mq halos to 5% for rich galaxy clusters of 
mass 10 15 /i _1 M in non-radiative simulations. Interestingly, the lower effective adiabatic index also increases the 
compressibility of the intrahalo medium, an effect that slightly increases the central concentration of the gas and 
the baryon fraction within the virial radius. We find that this can enhance the cooling rate onto central cluster 
galaxies, even though the galaxies in the cluster periphery become slightly less luminous as a result of cosmic ray 
feedback. 

Key words, galaxies: clusters: general - cosmic rays - methods: numerical. 



1. Introduction 

In recent years, the ACDM model has emerged as a 
highly successful 'concordance' model for cosmological 
structure formation. It conjectures that the dominant 
mass component in the universe consists of cold dark 
matter, and that a cosmological constant or dark en- 



ergy field add sufficient energy density to yield a spa- 
tially flat spacetime. This model is impressively success- 
ful in matching observational data on a large range of 
scales and epochs, includ ing the cosmic mic rowave back- 
ground fluc tuations fe.g.lSpergel et alll2003lh galaxy clus- 
tering (e.g. iTegmark et all l2004t ICole et all I2005Q . cos- 
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mic flows in the pres ent universe (e.g. IWillick et all Il997t 
iHudson et a~l 12004) . or the observational data on dist ant 
supernovae IjRiess et al 1 Il998t IPerlmutter et allll999j) . 

While the dynamics of the dark matter component in 
the ACDM model is now quite well understood and can 
be followed with h i gh accuracy in numerical simulations 
JPower et alll2003tlNa"varro et all l2004 iHeitmann et all 
120051) . the baryonic processes that regulate the formation 
of the luminous components of galaxies are much less well 
understood. Direct hydrodynamical simulations that fol- 
low the baryonic gas as well as the dark matter, face a 
number of 'small-scale' problems. For example, they tend 
to produce too many stars as a result of a 'cooling catas- 
trophe', unless effects like galactic outflows are included 
phenomenological way (e.g. ISpringel k, Hernauistl 
2003b|). The y lead to too concentrated disk galaxies (e.g. 



Abadi et al 1 EH) and fail to reproduce the observed 



shape of the luminosity function of gala xies in detail 
1 Murali et all l2002t iNagamine et all 12004) . 

By invoking strong feedback processes, semi-analytic 
models of galaxy formation are able to overcome 
these proble ms and to expl ai n a w id e array of galaxy 
properties llWhite fc FrenkL Il99li Kauffmann et al 



19931 iBaugh et all Il998t ISomerville fc PrimackL Il99£ : 



Cole et all l2000t ICroton et all l200fl . While this sup- 
ports the notion that feedback is crucial for the regulation 
of galaxy formation, it is unclear whether the physical 
nature of the feedback processes is correctly identified in 
the present semi-analytic models, or whether they merely 
give a more or less correct account of the consequences 
of this feedback. Direct hydrodynamic simulations can 
in principle be used to lift this ambiguity and to more 
directly constrain the physical processes at work. 

In most current models of galaxy formation, feed- 
back effects due to supernovae explosions and due to 
a photoionizing background are usually included, and 
more recently, some studies have consi dered quasar and 
radio activity by AGN as w ell (e.g. iDi Matteo et all 
200.4 ISiiacki fc Springe! boOfil) . However, perhaps surpris- 
ingly, magnetic fields and non-thermal pressure compo- 
nents from cosmic rays have received comparatively lit- 
tle attention thus far (with notable exceptions, includ- 
ing iKang etail Il996t iMiniatl 1200 it iMiniati etail l200ll 



mg 
iMin 



IMiniati l2002tlRvu & Kangj l2003ll2004j) . despite the fact 
that cosmic rays are known to contribute substantially to 
the pressure in the ISM of our own Galaxy. This is proba- 
bly at least in part due to the complexity of the cosmic ray 
dynamics, which when coupled to the galaxy formation 
process is very hard to describe analytically. Even when 
numerical methods are invoked, the cosmic ray physics is 
so involved that a number of simplifying approximations 
are required to make it tractable in a cosmological simu- 
lation setting, as we discuss here. 

In this study, our goal is to introduce the first cos- 
mological code of galaxy formation that treats cosmic 
rays self-consistently during the structure formation pro- 
cess. Our principal approach for capturing the cosmic 
ray physics has been laid out in a companion paper 



(En filin et all l2006|) . where we introduced a number of 
approximations to reduce the complexity of the problem. 
Fundamentally, we model the cosmic ray population of 
each fluid clement with a power law spectrum in particle 
momentum, characterized by an amplitude, a cut-off, and 
a fixed slope. Our model then accounts for adiabatic ad- 
vection of cosmic rays, and for injection and loss terms due 
to a variety of physical sources. Finally, we also include 
cosmic ray diffusion. The primary injection mechanisms 
we consider are supernova shocks and diffusive shock ac- 
celeration at structure formation shock waves. Since the 
efficiency of the latter is a sensitive function of the Mach 
number of the shock, we have also developed an on-the-fly 
shock finder for SPH calc ulations, which is descr ibed in a 
second companion paper ((Pfrommer et all l200fil) . 

In this paper, we use the theoretical model of 
Enfi lin et all l)2006() and cast it into a numerical formula- 
tion of cosmic ray physics that we im plement in the cosmo- 
logical TreeSP H code GADGET-2 (jspringel et all l200lt 
ISnringell . 120051) . We discuss our numerical approach in de- 
tail, including also various optimisations needed to keep 
the scheme efficient and robust. We then move on to show 
first results from applications of the model, ranging from 
isolated galaxies of different sizes, to cosmological simula- 
tions of galaxy cluster formation, and of homogeneously 
sampled boxes. Interestingly, cosmic rays can have a sub- 
stantial effect on dwarf galaxies, suppressing their star for- 
mation considerably. We show that this should leave a 
noticeable imprint in the luminosity function of galaxies, 
leading to a shallower faint-end slope. 

This paper is laid out as follows. In Section 2, we de- 
scribe the details of our implementation of cosmic ray 
physics, and in Section 3 we discuss our treatment of 
cosmic ray diffusion. Section 4 presents a number of test 
problems, which we used to verify the validity of results 
obtained by the code. We then describe in Section 5 a first 
set of simulations of isolated galaxies carried out with the 
new code. This establishes a number of principal effects 
found for the model. In Section 6, we then extend our 
analysis to more sophisticated, fully cosmological simula- 
tions of structure formation. We consider both galaxy clus- 
ters and dwarf galaxy formation at high redshift. Finally, 
Section 7 summarizes our conclusions and gives an out- 
look for future studies of cosmic rays in a cosmological 
context. 

2. Modelling cosmic ray physics 

In lEnfilin etail §006), we have introduce a new theoret- 
ical formalism for a simplified treatment of cosmic ray 
physics during cosmological structure formation. We also 
gave a detailed discussion of the physical background and 
the relative importance of various physical source and sink 
processes, and how they can be incorporated within the 
simplified framework. In this section of the present study, 
we describe the practical implementation of this model 
within the Lagrangian TreeSPH code GADGET-2, includ- 
ing also a concise summary of the those parts of the frame- 
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Fig. 1. Schematic illustration of the cosmic ray momen- 
tum spectrum in our two parameter model. We adopt a 
simple power-law description, where the slope of the cos- 
mic ray spectrum is given by a spectral index a, kept 
constant throughout our simulation. The normalization 
of the spectrum is given by the variable C, and the low- 
momentum cut-off is expressed in terms of a dimension- 
less variable q, in units of m p c where m p is the proton rest 
mass. 



work of lEnfilin et alJ l|200f)l) that we included in the code 
thus far. 

2.1. The cosmic ray spectrum and its adiabatic 
evolution 



As discussed in full detail in Enfilin et alJ l|2006l) . we as- 
sume that the cosmic ray population of each fluid element 
is made up of relativistic protons with an isotropic mo- 
mentum distribution function of the form 



f(p) 



dN 



Cp- a 6(p-q), 



(1) 



dpdV 

where C gives the normalization, q is a low momentum 
cut-off, and a is the power law slope. The momenta are 
expressed in dimensionless form in units of m p c, where 
m p is the proton mass. For the purposes of this paper, 
we will generally take a to be constant (a ~ 2.5 — 2.8), 
which should be a reasonable first order approximation in 
most cases relevant to galactic structure formation. We 
note however that a can in principle be made to vary in 
our formalism, at the price of a su bstantially increase d 
computational cost and complexity l(Enfilin et all Eflflfih . 
The pressure of this cosmic ray population is given by 



P C R = 



Ci 



■B 



2 3 



a 



(2) 



while the number density is simply tlcr = 
Cq l - a /(a-l).llem 



B n (a,b)= / x a - L (l-x) 
Jo 



6-1 



d.r 



(3) 
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Fig. 2. The function /3 a (q) introduced in equation (J5J, for 
several different values of the spectral slope a. 

denotes incomplete Beta functions. To describe the kinetic 
energy per cosmic ray particle for such a power-law pop- 
ulation we define the function 



Tenia, q) 



q a ~ 1 P a {q) + ^/l + l 



1 



which will be of later use. The quantity 

,' a — 2 3 — a 

P a {q) = B i 

1 + q 2 



(4) 



(5) 



is here introduced as a convenient abbreviation for the 
incomplete Beta function. We show (3 a (q) as a function of 
q for a few values of a in Figure El 

We here implement the cosmic ray model in a 
Lagrangian simulation code, where the advection of the 
cosmic ray population can be conveniently described sim- 
ply in terms of the motion of gas particles. In this ap- 
proach, the normalization of the spectrum should be ex- 
pressed in terms of a quantity normalized to mass, instead 
of the volume-normalized quantity C, with the translation 
between the two being simply given by the local gas den- 
sity p. In our case, it is convenient to absorb the proton 
mass into a redefinition of the amplitude, so we define 



C = C^ 



(6) 



as Lagrangian amplitude of the spectrum. Upon adiabatic 
changes of the gas density, the normalization of the spec- 
trum changes according to 



Co, 



while the momentum cut-off shifts as 



q(p) ( — I '/(i- 



(7) 



(8) 
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Here, we introduced a reference density po (for example 
set equal to the mean cosmic density) and a correspond- 
ing normalization Co and cut-off go a t this density. In 
our numerical implementation, we only have to follow the 
evolution of the adiabatic invariants Co and go due to 
non-trivial physical source and sink processes, releasing 
us from the task to compute adiabatic changes of the nor- 
malization and cut-off explicitly. Instead, they are simply 
accounted for by equations 10 and ©■ 

The number density ncR of relativistic CR protons can 
also be conveniently expressed in terms of the total proton 
density 1 , yielding 



n C R— = C— 
p a 



1 



= Cq- 



1 



(9) 



We can thus interpret n as something like a "cosmic ray 
to baryon ratio" . This quantity is an adiabatic invariant 
which can be followed accurately in dynamical simulations 
with our Lagrangian approach. 

Note that in our model we do not explicitely remove 
baryons from the reservoir of ordinary thermal matter 
when they are accelerated to become relativistic cosmic 
ray particles. This is a valid approximation, provided we 
have fi <C 1, which is always expected in our applications. 
In the calculations we performed so far, the fraction of 
protons contained in the relativistic phase typically re- 
mained far below a maximum value of n w 10~ 4 . The 
latter is already an exceptionally large value which we en- 
countered only in our most extreme tests, but it is still 
small enough that the reduction of the number density of 
thermal particles can be safely neglected. We note that 
cosmic ray confinement by magnetic fields holds in ideal 
magneto-hydrodynamics (MHD) when the mass fraction 
in relativistic particles is small. 

In a Lagrangian code, it is natural to express the cos- 
mic ray energy content in terms of energy per unit gas 
mass, e, which is given by 



C 



1 



1 



(10) 



Note that e refers to the energy normalized by the total 
gas mass, not by the mass of the cosmic ray particles alone. 
The specific energy content can also be expressed as 



T, 



CR ™CR 



CRn 



(11) 



In Figure |2| we show the distribution de/dlng of energy 
per logarithmic momentum interval, normalized to a spec- 
trum with vanishingly small cut-off. For spectral indices 
in the range 2 < a < 3, most of the energy is typically 
contained around q ~ 1, unless the cut-off of the actual 
spectrum lies higher than that, in which case the parti- 
cles just above the cut-off will dominate the total energy. 

1 We here loosly call p/m p the proton density. In our cosmo- 
logical applications, we of course use the mean particle mass 
where appropriate to account for the presence of heavier ele- 
ments and the ionization state of the gas. 
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Fig. 3. The distribution of cosmic ray energy per unit log- 
arithmic interval of proton momentum, for several differ- 
ent values of the spectral slope a. The distributions have 
been normalized to e(0) in each case. 

Due to our assumption that the momentum distribution 
extends as a power-law to infinity, the spectral index a is 
restricted to a > 2, otherwise the energy would diverge. 
For a < 3, the energy stays finite also for an arbitrarily 
low spectral cut-off. 

In our numerical scheme, every baryonic SPH parti- 
cle carries the adiabatic invariants go and Co as internal 
degrees of freedom for the description of the cosmic ray 
physics. These variables are then used to compute all phys- 
ical properties of the cosmic ray sector, as required for the 
force evaluations. For the gas dynamics, we are primarily 
interested in the effective pressure term due to the rela- 
tivistic particles, which are confined by the ambient mag- 
netic field. In our set of variables, this is compactly given 
as 

Cpc 2 



CR 



() 



-P a {q). 



(12) 



In the calculation of the hydrodynamic accelerations with 
the Euler equation, we can simply add this partial pres- 
sure due to CRs t o the ordinary thermal pressure (see 
Enfili n et al for further discussion) . This makes the 

interface with the ordinary hydrodynamical code conve- 
niently small, requiring only a small number of changes at 
well defined places. 

The effective adiabatic index of the cosmic ray pressure 
component upon local isentropic density changes is 



7CR 



d log Pgr 
d\ogp 



Zf3 a {q)Vl + q 2 



(13) 



On the other hand, when the pressure is expressed in terms 
of the cosmic ray energy density, we obtain 



Pgr 
pi 



(a-l)/3 a (q) 



3/3 Q (g) + 6g 1 -«(Vl + 9 2 -l) 



(14) 
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Fig. 4. Cosmic ray pressure in units of the cosmic ray en- 
ergy density, as a function of the spectral cut-off q. Except 
in the transition region from non-relativistic to relativistic 
behaviour, the cosmic ray pressure depends only weakly 
on q. In the ultra-relativistic regime, the ratio approaches 
Pcn/ipS) — (4/3 — 1), which is shown as the lower dot- 
ted line. The upper dotted line gives the expected value 
of (5/3 — 1) for an ideal gas. For the same energy content, 
cosmic rays always contribute less pressure than thermal 
gas. 



In Figure 21 we show the dependence of the right-hand- 
side of equation l|14|) on the spectral cut-off q, for differ- 
ent values of the slope a. For large values of q we obtain 
Pcn/ips) — (4/3 — 1), as expected for particles in the 
ultra-relativistic regime, while for low values of q the ratio 
is still significantly below the (5/3 — 1) expected for an 
ideal gas. However, it is clear that for given cosmic ray 
energy density, the pressure depends only weakly on the 
spectral cut-off; the value of e is hence much more impor- 
tant for the dynamics than the value of a. 

2.2. Including non-adiabatic CR processes 

The adiabatic behaviour of cosmic rays that are locally 
locked into the fluid by magnetic fields can be well traced 
with the above prescriptions. However, there are also a 
multitude of physical processes that affect the CR spec- 
trum of a gaseous mass-element in a non-adiabatic fash- 
ion. For instance, particles can be accelerated in strong 
shock waves to relativistic momenta and become cosmic 
rays. This process of diffusive shock acceleration should 
be particularly effective in high Mach number accretion 
shocks during cosmological structure formation, which can 
be traced by the hydrodynamical solver of our simulation 
code. On sub-resolution scales, violent shocks due to su- 
pernova explosions associated with stellar evolution may 
inject cosmic rays as well. Other astrophysical sources in- 



clude the ejection of high-energy particles in a jet from an 
accreting black hole. 

On the other hand, the cosmic ray population suffers 
a number of loss processes which will diminish the abun- 
dance over time if there is no new supply of freshly in- 
jected or accelerated protons. We shall here consider only 
the most prominent loss processes in the form of Coloumb 
losses that thermalize the cosmic ray energy, catastrophic 
losses that let the energy escape as radiation, and diffusion 
which washes out cosmic ray pressure gradients . 

As discussed above and in lEnfilin et all (2006). our cos- 
mic ray model requires three parameters to describe the 
state of the relativistic particle component of the gas. One 
parameter is the spectral index a, which we set to a con- 
stant value throughout the simulation volume, specified 
at the start of the simulation with a value motivated by 
the typically observed index in galactic systems. However, 
we do not restrict the range of non-adiabatic processes we 
consider to those with a similar injection index. Rather, we 
translate the injected cosmic ray properties into changes of 
amplitude and momentum cut-off within the framework of 
our simplified, fixed-slope model for the cosmic ray spec- 
trum. This translation is based on basic principles of mass 
and energy conservation. Despite this considerable simpli- 
fication, it is clear that the thermodynamic state of CR 
gas is considerably more complex than that of an ideal gas, 
where essentially everything is determined by the specific 
entropy alone. 

Given some changes within a simulation time step of 
the cosmic ray specific energy (de) and relative CR density 
(dn) associated with a fluid element, these changes can be 
cast into variations of the adiabatic invariants qo and Co 
of the cosmic ray population using a Jacobian matrix. We 
then obtain 



dC 



dC = C mpd f'I P / (g)d " (15) 



m p e - T p (q)n 



and 



dq = 



3 qo m p de-T CR dn 

Po J a - I m p e~ T p (q)n 



where we used the mean kinetic energy per cosmic ray 
particle, viz. Ten — sm p /n, as given by equation (0J, and 
defined 

„2 



T p (q) = Wl + q 2 -l)m p c 2 



(17) 



as the kinetic energy of a proton with normalized momen- 
tum q. Recall that Ten only depends on q and a, but does 
not depend on the normalization of the spectrum. 

However, this simple and fast translation scheme will 
only work for sufficiently small changes of the cosmic ray 
population. In our implementation, we therefore apply 
equations H15fl an d l|16fl only if the relative changes in 
cosmic ray energy and number density are less than a few 
percent. Otherwise, new cosmic ray spectral parameters 
in terms of Co and go are computed by explicitly solving 
equations and (|10(l . after applying the principles of 
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energy and particle conservation. While © can be easily 
solved for either q or Co, equation iJTUjl for the specific 
energy needs to be inverted analytically, but this can be 
done efficiently numerically ( s ee als o the discussion in sec- 
tion O and in lEnBlin et aill2006|) . 

Still, a naive application of energy and particle con- 
servation when adding a newly injected CR component to 
the current spectrum can cause unphysical results if the 
spectral cut-offs involved are very different. The reason lies 
in the strong dependence of the cosmic ray loss processes 
on particle momentum, together with our simplified spec- 
tral representation. As we will see, the life-time of cosmic 
ray particles grows monotonically with particle momen- 
tum. This dependence is particularly steep in the non- 
relativistic regime (ri osscs (p) ~ p 3 ), but becomes much 
shallower and eventually nearly flat in the mildly relativis- 
tic and strongly relativistic regimes. Simply injecting a CR 
component with very low cut-off to one with high cut-off 
while enforcing total energy and CR particle number con- 
servation will then result in a new composite spectrum 
where many of the original CR particles are represented 
as lower momentum particles. Consequently, their cool- 
ing times would be artificially reduced. Ultimately, this 
problem arises because the number density of the injected 
particles is dominated by low momenta, and these have 
cooling times much shorter than any relevant dynamical 
timescale. If this is the case, then it would make more 
sense to never inject this population to begin with, and 
to rather thermalize it instantly, thereby avoiding an un- 
physical distortion of the composite spectrum. 

The two injection processes we consider in this paper 
(shocks and supernova) both supply power-law distribu- 
tions of cosmic ray particles which start at very low ther- 
mal momenta. For them, we define an injection cut-off gi„j 
such that only the particles and the energy above the cut- 
off are injected, while the rest of the energy is instantly 
thermalized and added to the thermal reservoir. The cri- 
terion we choose for defining g; n j is 



^losses (^inj ) — ^"inj Qinj ) 

where Ti osscs (g) is the cooling timescale 



Tic 



|de/di| 



losses 



(18) 



(19) 



due to cosmic ray loss processes for a spectrum with 
cut-off at q and spectral slope a. As we will see, this 
timescale is independent of the normalization of the spec- 
trum, and inversely proportional to the density, provided 
the weak density-dependence of the Coulomb logarithm 
in the Coulomb loss-rate is neglected. Given the present 
cosmic ray energy content e, the timescale 



T inj {£, Qinj ) 



/ Q i„j(<7o,9inj) (de/di) 



(20) 



defines the current heating time due to the injection 
source, assuming that only the fraction /(go,ginj) above 
Qinj of the raw energy input rate (de/dt) i n j contributes effi- 
ciently to the build up of the cosmic ray population. Here 



go is the intrinsic injection cut-off of the source, which 
typically lies at very small thermal momenta, and a- m ^ is 
the slope of the source process. The factor f ai (g , gi„j) is 
given by 



/a inj (<7o,<7inj) = f~-J 



rCR.(«inj, <7inj) 

?CR(amj,go) 



(21) 



for ginj > go j otherwise by unity. 

Equation l|I8|) simply says that we only inject cos- 
mic ray particles at momenta where a spectral component 
could grow, given the rate of energy injection. It would 
be unphysical to assume that a spectrum develops that 
extends as a power-law to momenta lower than gj n j. The 
addition of this injection rule is hence necessary to make 
our simple model with a fixed spectral shape physically 
well-behaved. 

We note that equation l|I8|) will typically have two so- 
lutions, or may not have a solution at all if the current 
energy s in cosmic rays is high enough. In the former case, 
the physical solution is the smaller of the two, which lies at 
ginj < <Zmax, where g max is the place where the expression 
'Hosscs (9inj)/(go,9inj) attains its 



maximum, i.e. 



-j- [T lossos (g)/(g ,g)] 



9=<?n 



(22) 



which comes 



In the latter case, we choose q m j = g, 
closest to solving equation (|I8|) and naturally corresponds 
to the point where one expects the largest amount of cos- 
mic ray energy that can be present as a power-law with a 
balance between loss and source processes. 



2.3. Shock acceleration 

In our present model, we consider two primary sources for 
cosmic rays, diffusive shock acceleration and supernovae. 
In the former, a small fraction of the particles streaming 
through a shock front is assumed to diffuse back and forth 
over the shock interface, experiencing multiple accelera- 
tions This can only happen to particles in the high-energy 
tail of the energy distribution, and eventually results in a 
power-law momentum distribution function for the accel- 
erated particles. 

In the linear regime of CR acceleration, particles above 
a threshold momentum gi n j can be accelerated. They are 
redistributed into a power-law distribution in momentum 
that smoothly joins the Maxwell-Boltzmann distribution 
of the thermal post-shock gas. The slope of the injected 
CR spectrum is given by 

(23) 



Oil 



f 



where r 
shock ilBel 



/ p^ is t he density compression ratio at the 



fl978albh . If the shock is dominated by the 
thermal pressure, the spectral index can also be expressed 
through the Mach number M as 



4 - M 



37M 2 



2(M 2 - f ) 



(24) 
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The stronger the shock becomes, the flatter the spectrum 
of the accelerated CR particles, and hence the more high- 
energy particles are produced. Weak shocks on the other 
hand produce only rather steep spectra where most of the 
particles thermalize quickly. 

Due to the continuity between the power law and the 
thermal spectrum, the injected cosmic ray spectrum is 
completely specified by the injection threshold q , pro- 
vided the shock strength is known. We will assume that 
qo is at a fixed multitude x- m j of the average thermal post- 
shock momentum, p t h = \/2fcT2/ (m p c 2 ), i.e. qo = x ul] p t h- 
In this case, the fraction of particles that experience shock 
acceleration does not depend on the post-shock tempera- 
ture T2, and is given by 



Anii„ = 



'IT Cfe, 



(25) 



We will typically adopt a fixed value of x; n j ~ 3.5, mo- 
tivated by theoretic al studies of shocks in galactic super- 
nova remnants (e.g. iDrurv et all fl989'l. The fraction of 
injected supern oyae particles in st r ong shocks is then a 

r 4 Jd~ 
iKang fc .Tones! Il995h . 

In the linear regime of CR acceleration, the specific 
energy per unit gas mass in the injected cosmic ray pop- 
ulation is given by 



few times 10 4 llDrurv et all Il989t iJones fc KaneL 11993: 



Ae 



lin 



JcR(«inj, ginj) Ann 



(26) 



We can use this value to define a shock injection efficiency 
for CRs by relating the injected energy to the dissipated 
energy per unit mass at the shock front. The latter appears 
as extra thermal energy above the adiabatic compression 
at the shock and is given by Aiidiss = u 2 — Mir 7 , where 
u\ and U2 are the thermal energies per unit mass before 
and after the shock, respectively. The injection efficiency 
of linear theory is then given by 



A% 
Awdi, 



(27) 



However, the shock acceleration effect experiences satu- 
ration when the dynamical CR becomes comparable to 
the upstream ram pressure p\v\ of the flow. We account 
for thi s by adopting the limiter suggested by En filin et alJ 
(2006), and define as final acceleration efficiency 



1 — exp 



lin 



Cn 



(n 



(28) 



We will adopt C max = 0.5 for the results of this study. 
Thus, we take the injected energy to be 



Ae i: 



Cinj A^diss ? 



(29) 



and correspondingly, the injected particle number is given 
by 

An inj = T ^ A£inj (30) 

iCR.(ainj, go) 



where TcR(amj, <?o) is the mean kinetic energy of the ac- 
celerated cosmic ray particles. In practice, both Aej n j and 
An; n j will be lowered when we shift the actual injection 
point from q to q; n j , as determined by equation i|18|) , with 
the difference of the energies fed to the thermal reservoir 
directly. 

It is clear that the efficiency of CR particle accelera- 
tion depends strongly on the compression ratio, or equiva- 
lently on the Mach number of shocks. Interestingly, accre- 
tion shocks during cosmological structure formation can 
be particularly strong. Here we hence expect potentially 
interesting effects both for the forming intragroup and in- 
tracluster media, as well as for the intergalactic medium. 
However, in order to accurately account for the cosmic 
ray injection by structure formation shocks, we somehow 
need to be able to estimate the strength of shocks in SPH 
simulations. As SPH captures shocks with an artificial vis- 
cosity instead of an explicit shock detection scheme, this 
is a n on-trivial problem. 

In lPfrommer et all |2006), our second companion pa- 
per to this study, we have proposed a practical solution to 
this problem and developed a novel method for measuring 
the Mach number of shocks on the fly during cosmolog- 
ical SPH calculation. The method relies on t he en tropy 
formulation for SPH bv lSpringel fc Hernauis il J2002h . and 
uses the current rate of entropy injection due to viscosity, 
together with an approximation for the numerical shock 
transit time, to estimate the shock Mach number. The 
scheme works better than one may have expected, and it is 
in fact capable of producing quite accurate Mach number 
estimates, even for the case of composite gases with a ther- 
mal and a cosmic ray pressure component. In cosmologi- 
cal simulations, the method delivers Mach number statis- 
tics which agree well with results obtained with hydro- 
dynamical mesh codes that use explicit Riemann solvers 
i Rvn et a,ll|2003l:IPfrommer et alllioOfih . 



Having a reliable Mach number estimator solves an 
important problem when trying to account for cosmic ray 
injection by shocks in SPH. Another complication is posed 
by the shock broadening inherent in SPH, which implies a 
finite shock transit time for particles, i.e. a SPH particle 
may require several timesteps before it has passed through 
a shock and received the full dissipative heating. Note that 
the number of these steps depends on the timestep crite- 
rion employed, and can in principle be made very large 
for a sufficiently conservative choice of the Courant coeffi- 
cient. Unlike assumed in the above treatment of diffusive 
shock acceleration, we hence are not dealing with a dis- 
crete injection event, but rather need to formulate the cos- 
mic ray acceleration in a 'continuous fashion', in parallel 
to the thermal dissipation, such that the final result does 
not depend on how many timesteps are taken to resolve a 
broadened shock front. 

Fortunately, the above treatment can be easily ad- 
justed to these conditions. We can simply insert for At/diss 
in equation l|29l) the dissipated energy in the current 
timestep. This quantity is computed in the SPH formal- 
ism anyway, and in fact, we know that SPH will integrate 



8 



M. Jubelgas, V. Springel, T. Enfilin, C. Pfrommer: Cosmic rays in hydrodynamical simulations 



Aitdiss correctly through the shock profile, independent of 
the number of steps taken. This is because the correct pre- 
and post-shock state of the gas are enforced by the con- 
servation laws, which are fulfilled by the conservative SPH 
code. For the same reason, the Rankine-Hugoniot jump 
conditions are reproduced across the broadened shock. 
Note however that for computing the linear shock injection 
efficiency according to equation (J27J), we need to continue 
to use an estimate for the total energy dissipated across 
the shock, based on the Mach number finder. 

For simplified test calculation with the code, we have 
also implemented an option where we assume a constant 
injection efficiency of the shock acceleration process, along 
with a constant spectral index and momentum cut-off pa- 
rameter. Values for these parameters can then be cho- 
sen to represent the energetically most important types of 
shocks in the environment to be simulated. Such a simpli- 
fied injection mechanism can then also be used get an idea 
about the importance of the Mach-number dependence of 
the shock acceleration for different environments. 

2.4. Injection of cosmic rays by supernovae 

Strong shock waves associated with supernovae explosions 
are believed to be one of the most important cosmic ray 
injection mechanisms in the interstellar medium. However, 
similar to star formation itself, individual supernovae are 
far below our resolution limit in cosmological simulations 
where we need to represent whole galaxies, or more chal- 
lenging still, sizable parts of the observable universe. We 
therefore resort to a subresolution treatment for star for- 
mation and its regulation by supernovae, as proposed by 
ISpringel fe Hernauistl l)2003a(l . In this model, the interstel- 
lar medium is pictured as a multiphase medium composed 
of dense, cold clouds, embedded in a tenuous hot phase. 
The clouds form by thermal instability out of the diffuse 
medium, and are the sites of star formation. The massive 
stars of each formed stellar population are assumed to 
explode instantly, heating the hot phase, and evaporating 
some of the cold clouds. In this way, a tight self-regulation 
cycle for star formation in the ISM is established. 

To model the generation of cosmic rays, we assume 
that a certain fraction Csn — 0. 1 — 0.3 of the supernova en- 
ergy appears as a cosmic ray population l)Aharonian et all 
2006J iKang fc Jonesl 12006^1 . The total rate of energy in- 
jection by supernovae for a given star formation rate p+ 
depends on the IMF. Assuming a Salpeter IMF and that 
stars above a mass of 8 M Q explode as supernova with a 
canonical energy release of 10 51 ergs, we obtain roughly 
one supernova per 250 M Q of stellar mass formed, trans- 
lating to an energy injection rate per unit volume of csnP*, 
with e S N = 4 x 10 48 ergsMg 1 . We then model the CR en- 
ergy injection per timestep of a gas particle as 



AesN = Csncsn m* At, 



(31) 



where m* = p\j p is the particle's star formation rate per 
unit mass. Note that uncertainties in the IMF are not 



really important here as we have introduced a free param- 
eter, Csn, to control the amount of energy that is fed into 
cosmic rays. 

For the slope of the injected cosmic ray power-law we 
assume a p lausible value of cusn = 2.4 I Aharonian et all 
2004, 2006), and for the cut-off gsN, we can adopt the ther- 
mal momentum </sn = V '&2sn / (m p c 2 ) for a fiducial su- 
pernova temperature characteristic for the involved shock 
acceleration. Our choice of q;sn = 2.4 for the injection 
slope is motivated by the observed slope of ~ 2.75 in the 
ISM, and the realization that momentum dependent diffu- 
sion in a turbulent magnetic field with a Kolmogorov-type 
spectrum on small scales should steepen the injected spec- 
trum by p -1 / 3 in equilibrium. Our results do not depend 
on the particular choice for TgN, provided qsn -C 1. The 
change of the particle number density can then be com- 
puted with the mean kinetic energy Tcr(«sn, Qsn) of the 
injected power law. Using equations J3J and flllfl. this re- 
sults in 

A ~ CsNCSNm* 

An = m P \ A *- ( 32 ) 

J CR («SN , g SN) 

We not e that in the formalism of ISpringel fc Hernauistl 
(2003gj), we need to reduce the supernovae energy injected 
into thermal feedback (and to an optional wind model if 
used) by the fraction £sn that we assume powers cosmic 
ray acceleration. 

2.5. Coulomb cooling and catastrophic losses 

Charged particles moving through a plasma will gradu- 
ally dissipate their kinetic energy and transfer it to the 
surrounding ions and electrons by Coulomb interactions. 
The rate of this energy loss depends both on the physical 
properties of the surrounding medium and on the detailed 
momentum spectrum of the cosmic ray population. The 
latter in particular complicates an accurate determination 
of the Coulomb loss rate. 

Since particles with low momenta are most strongly 
affected by the Coulomb interactions, a qualitative conse- 
quence of this effect is that it induces a flattening of the 
spectrum; the high-momentum tail remains unchanged 
while the low-momentum cosmic ray particles dissipate 
their energy effectively to the thermal gas, and eventually 
drop out of the cosmic ray population altogether. 

In our model, we have deliberately abandoned a de- 
tailed representation of the cosmic ray spectrum of each 
fluid clement, in favour of the high computational speed 
and low memory consumption allowed by our simplified 
spectral model. We even opted to use a globally fixed spec- 
tral index, which means that we cannot represent a spec- 
tral flattening in detail. However, we can account for the 
effect of thermalization of the low-momentum particles in 
a simple and efficient way. To this end, we compute the 
energy loss by Coulomb-cooling over the entire spectrum, 
and then shift the low-momentum cutoff of our spectral 
model such that just the right amount of energy is removed 
in low-momentum particles, keeping the high-momentum 
part unaffected. 
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Fig. 5. The top panel shows the cooling times due to 
Coulomb losses (rising solid line) and hadronic dissipation 
(nearly horizontal line) as a function of the spectral cut- 
off. The dot-dashed line gives the total cooling time, while 
the vertical dashed lines marks the asymptotic equilibrium 
cut-off reached by the CR spectrum when no sources are 
present. The bottom panel shows the cooling time of ordi- 
nary thermal gas due to radiative cooling (for primordial 
metallicity) , as a function of temperature. The horizontal 
dashed line marks the cooling time of CRs with a high 
momentum cut-off (g > f), for comparison. In both pan- 
els, the times have been computed for a gas density of 
p = 2.386 x 10 -25 gcm -3 , which corresponds to the den- 
sity threshold for star formation that we usually adopt. 
Note however that the cooling times all scale as r oc 1/p, 
i.e. for different densities, only the vertical scale would 
change but the relative position of the lines would remain 
unaltered. 



More specifically, we follow lEnfilin et alJ |2006) and 
estimate the Coulomb loss rate of the CR population as 



fdt 
[dt. 



27rCe 4 



In 



2m c c 2 (/3p) 



Till; 



(33) 



f a - 2 



Here n c is the electron abundance, ui = yj A'Ke 2 n e /m e the 
plasma frequency, and (/3p) = 3Pcr/ (phc 2 ) is a mean 
value for our assumed spectrum. We also define a cool- 
ing timescale due to Coulomb cooling as 



|de/dt| ( 



(34) 



which depends only on the spectral cut-off q, and is in- 
versely proportional to density, modulo a very weak addi- 
tional logarithmic density-dependence though the plasma 
frequency. For a given energy loss Aec i n a timestep 
based on this rate, we can then estimate the corresponding 
change in cosmic ray number density as 



Anc = Aec 



T P (lY 



(35) 



i.e. we assume that the particles are removed at the low 
momentum cut-off q. From equations 115|) and (|I6[) . we can 
see that this will result in a gradual rise of the spectral cut- 
off q while the normalization will remain unchanged. The 
corresponding energy loss Aec is added to the thermal 
energy in the ordinary gas component, i.e. the Coulomb 
cooling process leaves the total energy content of the gas 
unaffected. We note that for large Coulomb cooling rates, 
we use an implicit solver to compute the new position of 
the spectral cut-off, leaving the amplitude parameter ex- 
actly invariant. This ensures stability even if the cut-off 
increases substantially in one timestep. 

Another class of loss processes for cosmic rays results 
occurs through inelastic collisions with atoms of the am- 
bient gas, resulting in the hadronic production of pions, 
which subsequently decay into 7-rays, secondary electrons, 
and neutrinos. In this case, the energy is ultimately dissi- 
pated into photons which tend to escape. So here the net 
effect is a loss of energy from the system, unlike in the 
case of Coulomb losses, where the dissipated cosmic ray 
energy heats the thermal reservoir. 

However, these 'catastrophic losses' can only proceed 
efficiently when the cosmic ray particles exceed the en- 
ergy threshold of gthrWpC 2 = 0.78 GeV for pion production 
(qt.hr = 0.83). The to tal loss rate can then be described by 
(Enfilin et~at l 120061) 



had 



cpa pp C ql a r C R(a,g jf ) 
2(a - I) m 2 



(36) 



where cf pp ~ 32 mbarn is the averaged pion production 
cross section, and g+ denotes q* = max {q, qthr}- The num- 
ber density of cosmic ray particles stays constant, however, 
due to conservation of baryon number in strong and elec- 
troweak interactions, i.e. we have Afihad = 0. This con- 
dition in turn implies that the changes of amplitude and 
cut-off of our spectral model due to this cooling process 
are related by 



AC 



Aq 
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is reached at the solution of the equation 



Fig. 6. Pressure of the cosmic ray population predicted 
for equilibrium between supernova injection on one hand, 
and Coulomb cooling and catastrophic losses on the other 
hand. The sold lines mark the pressure as a function of 
overdensity for two values of the injection efficiency £sn- 
The assumed threshold density for star formation, p — 



2.4 x 10 



" 25 _g_cm 3 



is derived from the multiphase model 



of ISprineel &: Hernauistl 1 2003aj) . The latter also predicts 
an effective equation of state for the star forming phase, 
shown as a dot-dashed line. The part below the threshold 
is an isothermal equation of state with temperature 10 4 K. 



As before, we also define a cooling timescale due to catas- 
trophic losses, which is given by 



Thad(<?) 



2m p T C R(a,q) 



\de/dt\ 



had 



cpo-pp Tenia, q*) 



(38) 

Note that Th a d becomes constant for q > qthi- 

In the top panel of Figure 03 we show the cooling 
timescales for Coulomb and catastrophic losses as a func- 
tion of the cut-off parameter q, at a fiducial density 
corresponding to our star formation threshold, assum- 
ing a spectral index a = 2.5. As expected, catastrophic 
losses dominate for high momentum cut-offs, and there- 
fore limit the lifetime of any cosmic ray population to 
r ~ 2 m p /(c/9(7pp), unless an injection process provides a 
resupply. For small cut-offs, Coulomb losses dominate and 
will rapidly thermalize the low momentum cosmic rays. 
The dot-dashed line shows the total loss timescale, defined 
by 

' ' (39) 



no 



3 (q) rciq) Thad(g) 



This timescale is monotonically rising with q. 

In the absence of any injection, the cosmic ray popula- 
tion will always evolve towards a fixed cut-off qa x , driven 
by the tendency of Coulomb cooling to increase the cut- 
off, while catastrophic losses tend to lower it. A balance 



1 



Tc(gfix) _ r C R(a,9fix) 



Thad ((/fix) Tp(<7fix) 



(40) 



which follows from equations (|15fl and Ijlfcifl . The vertical 
dotted line in Figure [5] marks this equilibrium point at 
(7fi x = 1.685 for a = 2.5. Once this fix-point is reached, 
only the spectral amplitude decays due to the cosmic ray 
cooling and dissipation processes. Finally, note that both 
the Coulomb cooling and the hadronic cooling time are 
inversely proportional to density. This also means that 
the fix-point qa x is density independent, but whether it 
can be reached in the available time at low density is a 
separate question. 

It is also interesting to compare the cosmic ray loss 
timescale to the thermal cooling timescale of primordial 
gas. The latter is also inversely proportional to density, 
but has in addition a strong temperature dependence. In 
the bottom panel of Figure 03 we show the thermal cool- 
ing timescale as a function of temperature, at the same 
fiducial density used in the top panel. For comparison, 
we show the asymptotic cosmic ray dissipation time scale 
(dashed line) , which is reached if the cosmic ray popula- 
tion is dominated by the relativistic regime. In this case, 
cosmic rays decay much slower than the thermal gas pres- 
sure in the intermediate temperature regime. This could 
be interesting for cooling flows in halos or clusters. Even 
a moderate cosmic ray pressure contribution in the dif- 
fuse gas in a halo of temperature T v „ ~ 10 5 — 10 7 K 
should tend to survive longer than the thermal gas pres- 
sure, which could influence the cooling rate. We will exam- 
ine this question explicitly in our numerical simulations of 
isolated halos. 

2.6. Equilibrium between source and sink terms 

The above considerations lead to an interesting question: 
What will the cosmic ray spectrum look like for a fluid el- 
ement at a density p high enough to allow star formation, 
such that is fed at a constant rate by supernova injection? 
We expect that after some time, a balance will be estab- 
lished between the supernova input on one hand and the 
cosmic ray losses due to Coulomb cooling and hadronic 
processes on the other hand. The energy content in the 
cosmic rays at this equilibrium point will then determine 
the CR pressure, and comparison of this pressure with the 
thermal pressure of the ISM will show whether 'cosmic ray 
feedback' could be important in regulating star formation 
in galaxies. 

To derive this equilibrium point, we first note that from 
the conditions Aesn + Aec + Aehad = and AnsN + 
Anc = 2 it directly follows that the mean energy per 
injected particle due to supernovae in equilibrium must 
be 



TcR(asN,9inj) = T p (q eq ) 



1 



Tc(<7eq) 



Thad(<Zcq) 



(41) 



Note that we always have Afihad = 0. 
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This is a relation between the effective injection cut-off of 
the supernova feeding, and the equilibrium cut-off q eq of 
the CR spectrum. In equilibrium, we also know that the 
supernova input will just balance the cooling losses for the 
spectrum with equilibrium cut-off q cq , yielding 

£ = n osscs (q cq ) f(q SN ,q hl j) . (42) 

On the other hand, our injection criterion of equation l|18|) 
tells us that e = ri osscs (g inj ) /(q SN , g inj ) (gf) SN , provided 
a solution for equation (|T%)| actually exists when the sys- 
tem is in the dynamic equilibrium. Assuming this for the 
moment, it follows that the CR loss timescales at q- ln j and 
q eq are equal, which in turn implies q- m j = q cq . In other 
words, the injection cut-off coincides with the cut-off of 
the equilibrium spectrum. The location of the equilibrium 
cut-off itself is given as solution of 

l + T C (g cq ) = ?CR.( a SN,gcq) ^ 
Thad(<Zeq) T p (q cq ) 

This is almost identical to equation l|4U|) which describes 
the fix-point of the cut-off if there is no injection. The 
only difference is the occurrence of cxsn instead of just a 
in the argument of the Tqb. function. The result of this 
will be a slight shift of the equilibrium position once the 
assumed spectral indices for supernova injection and the 
general cosmic ray population differ, while for q;sn = 
there will be no difference. At first sight it may seem sur- 
prising that the equilibrium position can be shifted away 
from the fix-point by an arbitrarily small supernova in- 
jection rate. However, recall that in the injection case we 
have described a truly invariant spectrum with a fixed am- 
plitude (which may take very long time to be established), 
while in the case without injection, the amplitude keeps 
falling on the cooling timescale. 

The above assumed that the injection condition (|18fl 
has a solution in equilibrium. This will be the case if 
q eq determined by equation 14311 is smaller or equal than 
9max as given by equation (|22l) . Otherwise the injection 
cut-off is given by (fr n j = q max , and the position of the 
equilibrium cut off is determined by equation i|42|) ■ A de- 
tailed comparison of the steady-state CR spectrum with 
and without our ap proximate description is provided by 
Enfi lin et alJ l|2006|) . There it is shown that dynamically 
important quantities like cosmic ray pressure and energy 
are calculated with ~ 10% accuracy by our formalism. 

Of primary importance for us is the equilibrium value 
of the cosmic ray energy content, because this will directly 
determine the CR pressure and hence the strength of po- 
tential feedback effects on star formation. The total cosmic 
ray energy injection rate by supernovae is related to the 
star formation rate by equation (|31|l . Once the equilibrium 
cut-off is q cq is known, the energy content in equilibrium 
will be given by equation l|42|) , so that the pressure is fully 
specified. For example, in the case of q cq < (? maX 7 the final 



pressure is therefore given by 

_ (a - l)^ a (<7eq) ^loss (<?eq) /q S n(9SN, qeq) r 

-tcr = ; <,sn£snP*- 

3/3 a (9eq) +6 (?e 1 q a ( V /l + 9e 2 q -l) 

(44) 

In our simulations, we combine the cosmic ray formal- 
ism with the subresolution multiphase model for the regu - 
lation of star formation bv lSpringel fc Hernauistl l)2003a|) . 
In the latter, the mean star formation rate is determined 
by the local density, and scales approximately as oc p 
above a density threshold for the onset of star formation. 
A detailed discussion of the multiphase model together 
with the precise de nsity dependence of the star forma- 
tion rate is given in ISpringel k Hernauistl j2003a|) . These 
authors also derive an effective equation of state for the 
ISM, which governs the assumed two-phase structure of 
the ISM. 

In Figure we show this effective pressure as a func- 
tion of density, using the parameters for gas consump- 
tion timescale, initial mas s function, and cloud evapora - 
tion efficiency discussed bv lSpringel k, Hernauistl l|2003a|) . 
which result in a star formation threshold of pth = 2.4 x 
lCT 25 gcm~ 3 . We also plot the expected cosmic ray pres- 
sure in the same diagram, for two different injection effi- 
ciencies of Csn = 0.1 and Csn = 0.3. Quite interestingly, 
the cosmic ray pressure exceeds the thermal pressure at 
the threshold density for star formation, but due to the 
shallow dependence of the equilibrium pressure on density, 
with approximately Pgr ~ P°' 5 (note that ri osse s oc p _1 ), 
we find that the pressure of cosmic rays only plays a role 
for the low density part of the ISM model, while regions 
with very high specific star formation rates should be at 
most weakly affected. 

For small efficiencies Csn < 0.01, we essentially expect 
no significant influence of cosmic rays from supernova on 
the regulation of star formation whatsoever. Even for the 
fiducial choice of Csn — 1, the influence would vanish for 
densities p > lOOpth- Based on these findings, we expect 
that galaxies which form most of their stars at compara- 
tively low densities should be potentially strongly affected 
by the cosmic ray feedback, while this influence should be 
weak or absent for vigorously star-forming galaxies with 
higher typical ISM densities. Our numerical results pre- 
sented later will confirm this picture. The fact that CR- 
pressure seems to become dominant around the star for- 
mation threshold may even suggest that cosmic rays could 
play an active role in defining this threshold. 

3. Cosmic ray diffusion 

Our treatment is based on the notion that the cosmic ray 
population is approximately locked into a gas fluid element 
by a locally present magnetic field. Even a weak ambient 
field makes the charged particle gyrate around the field 
lines, preventing them from freely travelling over macro- 
scopic distances with their intrinsic velocity close to to the 
speed of light. The cosmic ray particles may move slowly 
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along the field lines, but their perpendicular transport is 
strongly suppressed. 

However, occasional scattering of particles on mag- 
netic irregularities of MHD waves can displace the gy- 
rocentres of particles, such that a particle effectively 
"changes its field line" , allowing it to move perpendic- 
ular to the field. Since realistic magnetic field configu- 
rations are often highly tangled, or even chaotic, this 
can lead to sizable cross-field speeds of the particles. 
From a macroscopic point of view, the motion of the 
cosmic ray population can be described as a diffusion 
process, which is anisotropic with respect to the local 
magnetic field configuration. The theory of the respec- 
tive diffusion coefficients is complicated, and uncertain 
for t he case of turbulent magnetic field configu r ations 
e.g. iRechester fc Rosenbluthl Il978l iDuffv et all I l995| 
Bieber fe Matthaeusl Il997t iGiacalone fc Jokipii Il999t 
Naravan fc MedvedevU200lUEnfilinLl2003() . 

In principle, one could try to simulate the mag- 
netic field in SPH and then treat the diffusion in 
an anisotropic fashion. While recent advances in mod- 
elling MHD with SPH ar e promising l|Dolag et all Il999t 
IPrice fc Monaghanl Hp04) , these techniques still face se- 
vere numerical challenges when applied to simulations 
with radiative cooling. We therefore defer such an ap- 
proach to future work and model the diffusion isotropi- 
cally. Also, since we have no direct information about the 
local magnetic field strength and field configuration, we 
will invoke a phenomenological approach to estimate the 
effective diffusion coefficient as a function of the local con- 
ditions of the gas. It would be easy however to refine the 
spatial dependence of the diffusion coefficient once a more 
detailed physical scenario becomes available. 

Assuming isotropicy, the diffusion in the CR distribu- 
tion function f(p, x, t) can be written as 



% = V(«V/). 



(45) 



The diffusion coefficient k itself depends on the particle 
momentum (more energetic particles diffuse faster), and 
on the local magnetic field configuration. For definiteness 
we assume that the underlying MHD is turbulent with a 
Kolmogorov power spectrum, in which case the momen- 
tum dependence of k is given by 



(46) 



where we assumed relativistic particle velocities with v ~ 
c. Integrating the diffusion equation over particle mo- 
menta then yields 



dn 
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1 
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VKV(g3n CR ), 



(47) 



where q is the spectral cut-off. Because more energetic 
particles diffuse faster, we expect that the diffusion speed 
for the cosmic ray energy density will be a bit higher than 
for the number density itself. To account for this effect, we 



first multiply the diffusion equation with T p (p), and 
then integrate over the particle momenta. This results in 
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(48) 



where Sqr = ps is the cosmic ray energy density per unit 
volume. The factor Tcn(a — |, q)/Tcn(a, q) is larger than 
unity and encodes the fact that the diffusion in energy 
de nsity proceeds fast er than in particle number density. 
In lEnfilin et al.ll)2006h . we also give more general formulae 
for different power-law dependences of the diffusivity, and 
provide a more accurate treatment where the reduction of 
the diffusion rate at sub-relativistic energies is accounted 
for. 

3.1. Modelling the difFusivity 

Due to the lack of direct local information about the mag- 
netic field strength in our present numerical models, we 
parameterize the dependence of the diffusion coefficient 
on local gas properties in terms of a fiducial power-law 
dependence on the local gas density and gas temperature. 
In particular, we make the ansatz 



k 



(49) 



for the diffusion constant, which is effectively a three pa- 
rameter model for the diffusivity, specified by the overall 
strength kq of the diffusion at a reference density and tem- 
perature, and by the power-law slopes n p and tit for the 
density and temperature dependence, respectively. 

While clearly a schematic simplification, this parame- 
terization is general enough to allow an analysis of a num- 
ber of interesting cases, including models where the typical 
magnetic field strength has a simple density dependence, 
which can be a reasonable first order approximation for 
some systems , for example for th e diffuse gas in cluster 
atmospheres l|Dolag et alll2004b(l . 

For definiteness, we now construct such a very simple 
model, which will be the fiducial choice for the results on 
diffusion presented in this study. In Kolmogorov-like MHD 
turbulen ce, the domina nt parallel diffusivity is expected to 
scale as ijEnfilinl Eifll 



K (X Ib 



(50) 



where Ib gives a characteristic length scale for the mag- 
netic field of strength B. We start out by assuming that 
the magnetic energy density is a fixed fraction of the ther- 
mal energy content, which corresponds to 



B oc p V2 T i/2 



(51) 



An appropriate length scale for Ib is more difficult to es- 
timate, as it will sensitively depend on the level of local 
MHD turbulence, and on the build up of the magnetic 
field due to structure formation processes. For simplicity, 
we here assume that the length scale is related to the local 
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Jeans scale, which may be appropriate for the conditions 
of a multiphase interstellar medium where local density 
fluctuations constantly form clouds of order a Jeans mass, 
which are then in part dispersed by supernova-driven tur- 
bulence. This then gives a scaling of the form 



b oc p 



-1/2 j.l/2 



(52) 



Combining equations (|51|) and (|52|) , we obtain a model 
for the conductivity in the form k oc T 1 / 6 /? -1 / 2 , i.e. tit = 
1/6 and n p = —1/2. We fix the overall st rength by alludin 
to m e asurements i n our own Galaxy l|Berezinskii et aT 
Il990t ISchlickeiseii [2002) , which estimate a diffusivity 
along the magnetic field lines in the interstellar medium 
of approximately 



3 x 10 



27 



cm 



10 



kpc 2 
Gy7 



(53) 



Adopting our typical temperature and density values of 
the ISM as reference values, we end up with the following 
model for the diffusivity 



1/6 



k = 10 



kpc 2 (_T_\ 

Gyr Vl0 4 Ky Vl0 6 M Q kpc 



P 



-1/2 



(54) 



We will use this parameterization in the simulations with 
diffusion analysed in the this study. It is clear however that 
this model needs to be interpreted with a lot of caution, 
as the real diffusivity is highly uncertain, and may vary 
widely between different parts of the Universe. Improving 
the physical understanding of the strength of the diffusion 
will therefore remain an important goal for cosmic ray 
physics in the future. 

3.2. Discretizing the diffusion equation 

We still have to discuss how we numerically implement 
diffusion in our Lagrangi an SPH code. We h ere follow 
a similar strategy as in lJubelgas et alJ l|2004f) . where a 
new treatment of thermal conduction in SPH was dis- 
cussed and applied to simulations of clusters of galaxies. 
In essence, the very same techniques that can be used to 
solved the heat diffusion equation can also be applied to 
the cosmic ray diffusion needed here. 

We first rewrite the diffusion equations into a 
Lagrangian form that is matched to the variables evolved 
in our simulation code. This results in 



and 



p— = VnVD n , 
at 



where we have defined the abbreviations 

Q 



D K 



and 



1 1/3 T CR (a- 3,9) 
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D n = T pq l/i n 
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(55) 
(56) 

(57) 
(58) 



respectively. Our method for representing these equations 
in SPH is based on a discretisation scheme for the Laplace 
operat or that avoids second order derivativ es of the SPH 
kernel llBrookshawi 1198.4 iMonaghanl Il992ft . which makes 
the scheme robust against particle disorder and numerical 
noise. This gives us an evolution equation in the form of 



dsi \ - nij 2nij(D £ j — D £t i) 



dt 



PiPj 



(59) 



and similarly for the cosmic ray number density. We here 
introduced a symmetrization of the diffusivities according 
to 



= 2- 



(60) 



T" Kj 

based on the suggestion bv IClearv fc MonagharJ l)l999|) 
Furthermore, we replaced one of the D £ terms in the pair- 
wise diffusion term by the kernel interpolant 



Pk 



(61) 



As lJubelgas et al.l ( 2004ft show, such a mixed formulation 
between intrinsic particle variables and SPH-smoothed in- 
terpolants substantially improves the numerical stability 
against small-scale particle noise. The smoothing process 
suppresses strong small-scale gradients, while long-range 
variations and their diffusive evolution remain unchanged. 
Since we use an explicit time integration scheme, this 
behaviour prevents stability problems due to the typical 
'overshooting' problem that otherwise may arise due to 
strong local gradients from local outliers. 

Nevertheless, we still need to impose a comparatively 
strict timestep criterion to ensure proper integration of 
the diffusion. For the thermal conduction problem, we em- 
ployed a simple criterion that limits the relative change 
of thermal energy of a particle within a single timestep. 
Although the diffusion studied here is in principle very 
similar to the conduction problem, an equivalent criterion 
is not a good choice, simply because unlike for thermal 
conduction, the relevant reservoir can often be empty. In 
fact, we typically start simulations from initial conditions 
where all the cosmic ray particle densities are identical to 
zero. 

However, a closer look at the Green's function 
G(x, t) = (47TKt)~ 3 / 2 exp[-x 2 /(4fti)] of the diffusion pro- 
cess shows that differences between two points with a 
distance of |x| are diffused away with a characteristic 
timescale x 2 /k. Using the mean interparticle separation 
of SPH particles for the distance, this suggests the defini- 
tion of a diffusion timescale in the form 



1 / m 

tarn — — — 
K \ p 



2/3 



(62) 



We use this to limit the maximum timestep for particles 
to be constrained by 



1 / m 

At < e f diff = e- — 



2/3 



(63) 
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where we used e — 0.1 for the simulations presented in 
this study. This has provided us with a numerically stable 
cosmic ray diffusion while at the same time are prevented 
from becoming impractically small. 

As an additional refinement to the implementation of 
diffu sion, we have im plemented the method proposed by 
Ijubeleas et alJ l)2004|) to obtain a manifestly conservative 
scheme for cosmic ray energy and particle number, even 
when individual and adaptive timesteps are used. To this 
end, we rewrite equation Q59JI as 

3 

where we have defined a pairwise exchange term of cosmic 
ray energy in the form 

—rr = — | ,o —XijViWij. (65 

dt ,-,/-, |xy| 2 

In each system step, we now determine the change of the 
cosmic ray energy of particle i according to 

1 d 

rmtei = - £ At,(% - 5*)-^. (66) 

jk 

The double sum on the right can be simply evaluated by 
the ordinary SPH sums over the active particles, provided 
that for each neighbour j found for a particle i one records 
a change of AtiEij/2 for i, and a change of — AtiEij/2 for 
the particle j. We then apply the accumulated changes of 
the cosmic ray energy (or particle number) to all particles 
at the end of the step, i.e. not only to the ones that are 
active on the current timestep. In this way, we arrive at a 
scheme that manifestly conserves total cosmic ray energy 
and number density for each diffusive step. 

Finally, we note that we have implemented a further 
refinement in order to cope with technical problems asso- 
ciated with the situation in which isolated CR-pressurizcd 
particles are embedded in a background of particles with 
zero CR pressure. In the neighbourhood of such an iso- 
lated particle, the smoothed cosmic ray energy field D e j 
will be non-zero for particles that have themselves no CR 
component (yet). This can then in turn lead to exchange 
terms I£y between particles which both have zero cosmic 
ray pressure, leading to unphysical negative values for the 
energy in one of them. We found that this can be avoided 
if we limit the value of the interpolant D e ^ to be no more 
than a factor \ ~ 2.0 larger than the value D £t i for the 
particle itself. With this change, we recovered numerical 
stability of the diffusion in this situation. 

4. Numerical details and tests 

The numerical framework for cosmic ray physics presented 
above allows for very complex dynamics that interacts 
in non-linear and rather non-trivial manner with differ- 
ent aspects of ordinary hydrodynamics, and in particu- 
lar with the physics of radiative cooling and star forma- 
tion included in GADGET-2 to describe galaxy formation. 



Pitiit 




100 200 300 400 500 

x 



Fig. 7. Shock-tube test for a gas with thermal and cosmic 
ray pressure components. The simulation is carried out in 
a three-dimensional periodic box which is longer in the x- 
direction than in the other two dimensions. Initially, the 
relative CR pressure is Xcr = -PcR./-fth = 2 in the left 
half-space (x < 250), while we assume pressure equilib- 
rium between CRs and thermal gas for x > 250. The evo- 
lution then produces a Mach number M. = 3 shock wave. 
The numerical result of the volume averaged hydrodynam- 
ical quantities (p(x)}, (P(x)), (v x (x)), and (M(x)) within 
bins with a spacing equal to the interparticle separation of 
the denser medium is shown in black and compared with 
the analytic result shown in colour. 

Together with the nearly complete absence of analytic so- 
lutions, this makes a direct validation of the numerical 
implementation of cosmic ray physics in the simulation 
code particularly hard. 

However, there are a few areas where the careful checks 
of individual subroutines of the code that we carried out 
can be augmented with tests of problems where analytical 
solutions are known. One such area are hydrodynamical 
shock waves that involve a cosmic ray pressure compo- 
nent. This allows us to test one of the most interesting 
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dynamical aspects of cosmic ray physics which is the in- 
troduction of a variable adiabatic index 7. Note that the 
pressure of a 'hybrid gas' with a thermal and a cosmic ray 
energy density can no longer be described with the simple 
parameterizations used for a polytropic gas. In particular, 
both the relative contribution of the cosmic ray pressure 
and the adiabatic index of the CR pressure component 
itself will change during an adiabatic compression. While 
more complex than for an ideal gas, the Riemann problem 
for a sho cks in such a composit e gas can be solved ana- 
lytically ijPfrommer et all l2006|) . and we will use this as 
a test for our numerical treatment. 

Another aspect which can be tested with simple toy 
set-ups is the diffusion of cosmic rays. To this end we will 
consider a geometrically simple initial cosmic ray distri- 
bution, together with the gas being forced to be at rest. 
This allows us to test the correct diffusion speed, and the 
conservative properties of the diffusion process. 

4.1. Implementation issues 

We note that the treatment of a new physical process in 
SPH often requires a modification of the timestep criterion 
to ensure proper time integration of the added physics. In 
the case of cosmic rays, it turns out that an additional 
limit on the timestep is not really required, because the 
Courant criterion for hydrodynamics is already general 
enough and automatically adjusts the timestep if the cos- 
mic ray dynamics requires it, as the latter then induces a 
large sound speed due to extra cosmic ray pressure. The 
non-adiabatic source and sink terms on the other hand are 
essentially encapsulated in subresolution models, and are 
therefore not particularly demanding with respect to the 
integration timestep. 

As can be seen in equations (|1UII and (|12|) . calculating 
physical properties like pressure and specific energy of the 
cosmic ray component involves the evaluation of incom- 
plete beta functions, which can be rather expensive nu- 
merically. Also note that in the individual timestep scheme 
of GADGET-2, which is essential for simulations with a 
large dynamic range, the pressure for all SPH particles 
needs to be set at every system step, even if a particle is 
'passive' and does not receive a force computation itself 
in the current step. Clearly, a costly evaluation of special 
functions for the pressure would therefore imply a signifi- 
cant burden in terms of processor time. We have therefore 
implemented a series of look-up tables discretized in log q 
that allow us a fast evaluation of terms involving costly 
incomplete beta functions. Interpolating from these ta- 
bles allows a rapid and accurate computation of all cosmic 
ray related quantities without special function evaluations 
during the simulation. 

The same technique is also used to numerically invert 
the specific energy Tor of equation (|10l) for q, a task that 
arises when the mean CR energy is updated. Here we again 
resort to a pre-computed look-up table in which we lo- 
cate the new mean CR energy, and then interpolate in the 



table to determine the spectral cutoff q. Once the spec- 
tral cutoff and injected CR-to-baryon fraction is found, 
the new spectral normalization C can be computed from 
equation @. As a final step, we can then update the ef- 
fective hydrodynamic pressure due to cosmic rays of the 
particle in question. 

For treating Coulomb losses in an accurate and ro- 
bust way, we use an implicit scheme because of the very 
sensitive dependence of the cooling rate on the spectral 
cut-off q. In fact, in order to ensure that this cooling pro- 
cess leaves the normalization of the spectrum constant and 
only increases the amplitude q, we solve the following im- 
plicit equation 

e(q',C) = e(q,C)- £ -^f^At (67) 

for a new spectral cut-off q' when the cooling lasts for a 
time interval At. This scheme is robust also in cases where 
At exceeds Tc(q). Unlike the Coulomb losses, the timescale 
of hadronic losses is comparatively long and varies little 
with the spectral cut-off, as seen in Figure 03 It therefore 
is easier to integrate accurately and does not require an 
implicit solver. 

4.2. Shocks in cosmic ray pressurized gas 

We performed a number of three-dimensional test sim- 
ulations that follow a shock wave in a rectangular slab 
of gas, which is extended along one spatial dimension. 
Periodic boundary conditions in the directions perpendic- 
ular to this axis are used to make sure that no boundary- 
effects occur. This allows us to simulate a planar shock 
in 3D which can then be compared to a corresponding 
one-dimensional analytic solution. The initial conditions 
of our shock-tube tests were set-up with relaxed 'glass' 
structures of particles initially at rest, and by giving the 
two halves of the slab different temperatures and cosmic 
ray pressures. The particle mass was constant and cho- 
sen such that the mean particle spacing was 1 length unit 
in the high-density regime. While a reduction of the par- 
ticle mass in the low-density region would have resulted 
in an increased spatial resolution there, our constant par- 
ticle mass set-up is more appropriate for the conditions 
encountered in cosmological simulations. 

In Figured we show the state of the system after a 
time of t = 0.3, for an initial density contrast of 5, a to- 
tal pressure ratio of 35.674, a homogeneous mixture of 1/3 
cosmic ray pressure and 2/3 thermal pressure contribution 
on the left-hand side, and pressure equilibrium between 
CRs and thermal gas on the right-hand side. A shock with 
Mach number M = 3 travelling to the right, a rarefaction 
wave to the left, and a contact discontinuity in between 
develop for these initial conditions. The analytically pre- 
dicted shock front position and density distribution are 
matched nicely by the simulation. Due to the smooth na- 
ture of SPH simulations, the density jump at the shock 
at x ~ 430 is not a sharp discontinuity, but stretches over 



16 



M. Jubelgas, V. Springel, T. Enfilin, C. Pfrommer: Cosmic rays in hydrodynamical simulations 



mm 



f= 0.0 Gyr _ 



f=0. 1 Gyt 





' 

1 — 0.2 Gyt 


:. ' 


i . 




i 





1 1 L 1 1 J 

f= 0.5 Gyr _ 


\ 


1 


V " 







' "■" 

f = 1.0 Gyr _ 


: \ 


■ • ' 


v ; 


. . i 






f = 2.0 Gyr _ 


- } 


' 


V " 



Fig. 8. Time evolution of a step function in cosmic ray energy density due to diffusion (the spectral cut-off and slope 
of the cosmic rays are constant throughout the volume). The population of SPH particles is kept at rest. The times 
shown in the different panels are (from top left to bottom right): t = 0, O.f, 0.2, 0.5, f.O and 2.0 Gyr. Black dots give 
particle values at the corresponding time, while the red line shows the analytical solution. The diffusivity k is constant 
at lkpc 2 Gyr -1 on the left hand side, and four times higher on the right hand side. 



a small number of mean interparticle spacings. The con- 
tact discontinuity at x ~ 375 is reproduced well, with 
only a small 'blib' seen in both the density and the pres- 
sure profile, which is characteristic for SPH in shock tube 
tests. Note that cosmic ray pressure dominates over ther- 
mal pressure on the left side of the contact discontinuity 
due to adiabatic rarefaction of the initially CR-dominated 
state on the left-hand side, while this is reversed on the 
other side because CRs are adiabatically compressed at 
the shock while the thermal gas experiences entropy in- 
jection. The rarefaction wave traveling to the left shows 
the expected behaviour over most of its extent, only in the 
leftmost parts at x ~ 100 some small differences between 
the analytic and numerical solution are seen. However, the 
overall agreement is very reassuring, despite the fact that 
here a shock in a composite gas was simulated. 

This shows that the simulation code is able to correctly 
follow rapid compressions and rarefactions in a gas with 
substantial cosmic ray pressure support, including shocks 
that feed their dissipated energy into the thermal com- 
ponent. In cosmological simulations where diffusive shock 
acceleration of cosmic rays is included, some of the dissi- 
pated energy will instead be fed into the cosmic ray reser- 
voir, so that there the resulting shock behaviours can be 
yet more comp licated. We also note th at our shock detec- 
tion technique l|Pfrommer et all f2006) is able to correctly 
identify the shock location on-the-fly during the simula- 
tion, and returns the right Mach number in the peak of the 
shock profile, where most of the energy is dissipated. We 



can use this to accurately describe the Mach-number de- 
pendent shock-injection efficiency of cosmic rays in shocks. 

4.3. Cosmic ray diffusion 

In order to test the diffusion part of the code in a clean 
way, we use a system of gas particles that are at rest, which 
avoids the complications that would otherwise occur due 
to the motions of particles. We achieve this by setting all 
particle accelerations to zero, such that the gas in fact 
behaves essentially like a solid body. As a side effect, the 
densities remain constant over time, such that all varia- 
tions in the distribution of cosmic rays are entirely due to 
diffusive transport (we also switched of the Coloumb and 
catastrophic losses for this test). For this idealized situa- 
tion, analytic solutions for the diffusion problem can be 
derived which can then be readily compared with numer- 
ical results. 

For definiteness, we set up a periodic slab of matter 
with density 1 x 10 10 M Q kpc -3 , spanning a basic volume 
of 10 x 10 x 100 kpc. The periodicity across the short di- 
mensions ensures the absence of boundary effects, such 
that we can compare the numerical results to effectively 
one-dimensional analytic solutions. The cosmic ray distri- 
bution was initialized such that the energy density due to 
relativistic particles has a sharp step. The spectral cutoff 
at both sides of the step was set equal to q = 0.3 in the test 
discussed here, but we note that similar results are also 
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obtained for different choices. Again, we used an irregu- 
lar glass-like configuration as initial particle distribution 
in order to more realistically model the noise properties 
in density fields encountered in cosmological applications. 
Note that small-scale numerical no ise can be prob l emati c 
for the treatment of diffusion fe.g. Ijubeleas et all l2004j) . 
so this is an important aspect for testing the robustness 
of the scheme. 

In real physical applications, the diffusion implemen- 
tation will have to deal with a spatially varying diffusion 
coefficient. In particular, there will be steep gradients in 
the diffusivity at phase transition between the cold, dense 
gas and the hot, yet thin ambient intergalactic and intra- 
cluster medium. It is therefore advisable to verify that the 
implemented numerical scheme for the diffusion is well be- 
haved at sharp jumps of the diffusivity. We incorporate 
this aspect into our test scenario by setting up a fiducial 
temperature-dependent cosmic ray diffusivity of 



f.O 



kpc 2 
Gy7 



1 000K j 



(68) 



and varying the gas temperature from 1000K in the left 
half of the matter slab to 4000K in the right, causing an 
increase of the diffusivity by a factor of 4 across the x = 
plane. Of course, this particular choice of conductivity 
is arbitrary, but the chosen values are not too dissimilar 
compared with what we will encounter in cosmological 
simulations later on. 

In Figure|Sl we present the time evolution for diffusion 
of an initial step function. The simulation was run over a 
time span of 2 Gyr, and for a number of times in between, 
we compare the spatial distribution of the cosmic ray en- 
ergy d ensity obtained numeric ally with the analytical so- 
lution ( Pfrommer et all l2006|) for the problem (shown in 
red). The match of the numerical result and the analytic 
solution is very good, especially at late times. In fact, after 
t = I Gyr, we no longer see any significant deviation be- 
tween the numerical solution and the analytical one. The 
code reliably traces the flattening of the cosmic energy 
density jump over time. The largest differences occur in 
the very early phases of the evolution, at around the initial 
discontinuity, as expected. Due to the smoothing inherent 
in SPH and our diffusion formulation, sharp gradients on 
very small-scales are washed out only with some delay, 
but these errors tend to not propagate to larger scales, 
such that the diffusion speed of large-scale gradients is 
approximately correct. 

Note that small-scale noise present in the initial cos- 
mic ray energy distribution is damped out with differ- 
ent speeds in the left and right parts of the slab. This is 
due to the different conductivities in the low- and high- 
energy regimes, which give rise to characteristic diffusion 
timescales of i<jiff = 1 Gyr and tag = 0.25 Gyr, respec- 
tively, for our mean interparticle separation of I kpc, con- 
sistent with Eqn. (|62(l . We note that we have verified the 
good accuracy of the diffusion results for a wide range of 
matter densities and diffusivities, including also cases with 



stronger spatial variations in diffusivity. We are hence con- 
fident that our numerical implementation should produce 
accurate and robust results in full cosmological simula- 
tions, where the diffusivity can show non-trivial spatial 
dependences. 

5. Simulations of isolated galaxies and halos 

We now turn to a discussion of the effects of our cosmic 
ray model on the galaxy formation process. Due to the 
complexity of the involved physics, which couples radia- 
tive cooling, star formation, supernova feedback, cosmic 
ray physics, self-gravity, and ordinary hydrodynamics, it 
is clear however that our analysis cannot be fully exhaus- 
tive in this methodology paper. Instead, our strategy is to 
provide a first exploration of the most important effects 
using a set of simulations with idealized initial conditions, 
and a restricted set of full cosmological simulations. This 
can then guide further in-depth studies of the individual 
effects. 

One of the possible effects of cosmic ray physics is that 
the injection of CRs due to supernovae may alter the reg- 
ulation of star formation by feedback, which may directly 
translate into observable differences in forming galaxies. 
Since CR-pressurized gas has a different equation of state 
than ordinary thermal gas, it may rise buoyantly from 
star-forming regions, which could perhaps help to produce 
outflows from galactic halos. Also, because energy stored 
in cosmic rays will be subject to different dissipative losses 
than thermal gas, we expect that the radiative cooling of 
galaxies could be altered. Of special importance is also 
whether the strength of any of these effects shows a de- 
pendence on halo mass, because a change of the efficiency 
of galaxy formation as a function of halo mass is expected 
to modify the shape of the resulting galaxy luminosity 
function. 

Another intriguing possibility is that the total bary- 
onic fraction ending up in galactic halos could be modi- 
fied by the additional pressure component provided by the 
relativistic particle population. In particular, the softer 
equation of state of the cosmic-ray gas component (its 
adiabatic index varies in the range 4/3 < 7cr < 5/3) 
could increase the concentration of baryonic matter in 
dark matter potential wells, because the pressure increases 
less strongly when the composite CR/thermal gas is com- 
pressed. On the other hand, a partial cosmic ray pressure 
support might reduce the overall cooling efficiency of gas 
in halos, causing a reduction of the condensated phase of 
cold gas in the centres. 

To examine the non-linear interplay of all these effects, 
we will study them in a number of different scenarios. We 
first use isolated galaxy models which allow a precise con- 
trol over the initial conditions and an easy analysis and 
interpretation of the results. Next, we use non-radiative 
cosmological simulations to investigate the efficiency of 
CR production at structure formation shock waves. We 
then use high-resolution cosmological simulations that in- 
clude radiative cooling and star formation to study the 
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Fig. 9. Time evolution of the star formation rate in isolated halos of different mass which are initially in virial 
equilibrium. In each panel, we compare the star formation rate in simulations without cosmic ray physics (solid red 
line) to two runs with different injection efficiency of cosmic rays by supernovae, Csn = 0.1 (blue lines) and £sn = 0.3 
(green lines), respectively. From top left to bottom right, results for halos of virial mass 1O 9 /i _1 M0 to 10 12 h~ 1 M & 
are shown, as indicated in the panels. Efficient production of cosmic rays can significantly reduce the star formation 
rate in very small galaxies, but it has no effect in massive systems. 



formation of dwarf galaxies, with the aim to see whether 
our identified mass trends are also present in the full cos- 
mological setting. We also use these simulations to inves- 
tigate whether CRs influence the absorption properties of 
the intergalactic medium at high redshift. Finally, we use 
high-resolution 'zoomed' simulations of the formation of 
clusters of galaxies to study how CR injection by accre- 
tion shocks and supernovae modifies the thermodynamic 
properties of the gas within halos. 



5.1. Formation of disk galaxies in isolation 



As a simple model for the effects of cosmic ray feedback 
on disk galaxy formation, we consider the time evolution 
of the gas atmospheres inside isolated dark matter halos. 
The initial conditions consist of a dark matter potential 
with a structure motivated from cosmological simulations, 
combined with a hydrostatic gas distribution initially in 
equilibrium within the halo. We then consider the evolu- 
tion of this system under radiative cooling, star formation 
and cosmic ray production by supernovae. We expect that 
the gas in the centre looses its pressure support by cool- 
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Fig. 10. Phase-space diagram of the star-forming phase in two simulations with halos of different mass. In these 
fiducial simulations, we included cosmic ray physics but ignored the cosmic ray pressure in the equations of motion, 
i.e. there is no dynamical feedback by cosmic rays. However, a comparison of the cosmic ray pressure and the thermal 
pressure allows us to clearly identify regions where the cosmic rays should have had an effect. For graphical clarity, 
we plot the pressures in terms of a corresponding effective temperature, T c s = (fi/k)P/p. Above the star formation 
threshold, the small galaxy of mass 10 9 h Mq shown in the left panel has a lot of gas in the low-density arm of the 
effective equation of state, shown by the curved dashed line. On the other hand, the massive 10 12 h~ 1 M Q galaxy shown 
on the right has characteristically higher densities in the ISM. As a result, the cosmic ray pressure is insufficient to 
affect this galaxy significantly. Note that the falling dashed line marks the expected location where cosmic ray loss 
processes balance the production of cosmic rays by supernovae. We show the systems at time t — 2.0 Gyr after the 
start of the evolution. 



ing, and then collapse s into a rotationally supp orted disk 
that forms inside-out IjFall fc EfstathiouLll980|) . 

It is clear that this is a highly idealized model for disk 
galaxy formation, which glosses over the fact that in a 
more realistic cosmological setting galaxies originate in a 
hierarchical process from the gravitational amplification 
of density fluctuation in the primordial mass distribution, 
gradually growing by accretion and merging with other 
halos into larger objects. However, the simplified approach 
we adopt here should still be able to capture some of the 
basic processes affecting this hierarchy, and it does so in 
a particular clean way that should allow us to identify 
trends due with galaxy mass due to cosmic rays. 

We model the dark matter and baryonic con- 
tent of our isolated halos a s NFW density profiles 
ijNavarro. Frenk fc White! . Il996|) , which we slightly soften 
at the centre to introduce a core into the gas density, with 
a maximum density value lying below the threshold for 
star formation, allowing for a 'quiet' start of the simula- 
tions. The velocity dispersion of the dark matter and the 
temperature of the gas were chosen such that the halos 
are in equilibrium initially, i.e. when evolved without ra- 
diative cooling, the model halos are perfectly stable for 
times of order the Hubble time. We also impart angular 
momentum onto the halo with a distribution inside the 



halo consistent with results obtained from full cosmologi- 
cal simulations (|Bullock et all 12001 1. 

We simulated a series of host halos with masses varying 
systematically between 10 9 /i _1 M Q and 10 12 ft _1 M . In all 
cases, we adopted a baryon fraction of fib/f^m = 0.133, 
and a matter density of f2 m = 0.3. We typically repre- 
sented the gas with 10 5 particles and the dark matter with 
twice as many. In some of our simulations, we also replaced 
the live dark halo with an equivalent static dark matter 
potential to speed up the calculations. In this case, the 
contraction of the dark matter due to baryonic infall is 
not accounted for, but this has a negligible influence on 
our results. We have kept the concentration of the NFW 
halos fixed at a value of c — 12 along the mass sequence, 
such that the initial conditions are scaled versions of each 
other which would evolve in a self-similar way if only grav- 
ity and ideal hydrodynamics were considered. However, 
this scale-invariance is broken by the physics of cooling, 
star formation and cosmic rays. 

When one of these halos is evolved forward in time, 
radiative cooling leads to a pressure loss of the gas in the 
centre, which then collapses and settles into a rotationally 
supported cold disk. In the disk, the gas is compressed 
by self-gravity to such high densities that star formation 
ensues. Unfortunately, the physics of star formation is 
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not understood in detail yet, and we also lack the huge 
dynamic range that would be necessary do directly fol- 
low the formation and fragmentation of individual star- 
forming molecular clouds in simulations of whole galax- 
ies. In this study, we therefore invoke a sub-resolution 
treatment for star fo r mation , in the form described by 
ISpringel k, Hernauistl l)2003a|K The model assumes that 
the dense interstellar medium can be approximately de- 
scribed as a two-phase medium where cold clouds form 
by thermal instability out of a diffuse gaseous phase. The 
clouds are the sites of star formation, while the super- 
novae that accompany the star formation heat the diffuse 
medium, and, in particular, evaporate some of the cold 
clouds. In this way, a self-regulation cycle for star forma- 
tion is established. 

When our new cosmic ray model is included in our sim- 
ulation code, a fraction of the deposited supernova energy 
is invested into the acceleration of relativistic protons, and 
hence is lost to the ordinary feedback cycle. While this en- 
ergy no longer directly influences the star formation rate, 
it has an indirect effect on the star-forming gas by provid- 
ing a pressure component that is not subject to the usual 
radiative cooling. If this pressure component prevails suf- 
ficiently long, it can cause the gas to expand and to lower 
its density, thereby reducing the rate of star formation. 

In FigureE! we show the time evolution of the star for- 
mation rate for four different halos masses, ranging from 
10 9 ft. _1 M to 10 12 h~ 1 M Q . For each halo mass, we com- 
pare three different reference simulation where the 
ordinary model of ISnringel fc HernauistJ (jifl03a) without 
cosmic rays was used, and two simulations where cosmic 
ray production by supernovae was included (without al- 
lowing for diffusion), differing only in the assumed effi- 
ciency of Csn = 0.1 and £sn = 0.3 for this process, respec- 
tively. Interestingly, the simulations with cosmic rays show 
a substantial reduction of the star formation rate in the 
two small mass systems, but already for the 10 11 /i~ 1 M Q 
halo the effects becomes comparatively small, while for 
the massive halo of mass 10 12 /i _1 M0, no significant dif- 
ferences can be detected. Clearly, the ability of cosmic 
ray feedback to counteract star formation shows a rather 
strong mass dependence, with small systems being af- 
fected most. For higher efficiencies £sn of CR-production 
by supernovae, the reduction of the star formation rate 
becomes larger, as expected. 

Figure ITUI provides an explanation for this result, and 
also elucidates the origin of the oscillatory behaviour of 
the SFR in the CR-suppressed cases. In the figure, we 
show phase-space diagrams of the gas particles of the 
10 9 /i _1 M and 10 12 /t _1 M halos, respectively, in a plane 
of effective temperature versus density. We plot the ther- 
mal pressure and the cosmic ray pressure separately. In 
order to cleanly show whether a dynamical effect of cos- 
mic rays can be expected, we here use a fiducial simula- 
tion where the cosmic ray pressure is ignored in the equa- 
tions of motion, but is otherwise computed with the full 
dynamical model. As Figure EH demonstrates, the bulk 
of the star-forming gas in the massive halo lies at much 



higher density and higher effective pressure than in the 
low mass halo. Because the cosmic ray pressure exceeds 
the effective thermal pressure of the multi-phase ISM only 
for moderate overdensities relative to the star formation 
threshold, most of the gas in the 10 12 /i -1 Mq halo is sim- 
ply too dense to be affected by the cosmic ray pressure. 
We note that the relative sizes of the two pressure compo- 
nents are consistent with the analytic expectations shown 
in Figure El In fact, these expectations are replicated as 
dashed lines in Figure an d are traced well by the bulk 
of the particles. Because the shallower potential wells in 
low-mass halos cannot compress the gas against the effec- 
tive pressure of the ISM to comparably high overdensities 
as in high-mass halos, it is therefore not surprising that 
the cosmic ray pressure becomes dynamically important 
only in small systems. 

Figure ITUI also makes it clear that in the regime where 
cosmic ray pressure may dominate we cannot expect a dy- 
namically stable quasi-equilibrium with a quiescent evolu- 
tion of the star formation rate. This is simply due to the 
decline of the effective cosmic ray pressure with increas- 
ing density of the ISM, a situation which cannot result 
in a stable equilibrium configuration where self-gravity is 
balanced by the cosmic ray pressure. Instead, the system 
should be intrinsically instable in this regime. When some 
gas becomes dense enough to start star formation, it will 
first have no cosmic ray pressure support but it will be 
stabilized against collapse by the thermal pressure of the 
ISM that is quickly established by supernova feedback. 
After some time, the ongoing star formation builds up a 
cosmic ray pressure component, which eventually starts to 
dominate, at which point the gas is driven to lower den- 
sity. As a result, the star formation rate declines strongly. 
After some time, the CR pressure is dissipated such that 
the gas can collapse again. Star formation will then start 
again and the 'cycle' can repeat. This scenario schemat- 
ically describes the origin of the oscillations in the star 
formation rate seen in the results for the 10 9 /i _1 M Q and 
10 10 /i _1 M Q halos when cosmic rays are included. 

Another view of the halo mass dependence of the ef- 
fects of cosmic ray feedback on star formation is given 
by Figure El Here we show the integrated stellar mass 
formed up to time t — 3 Gyr, normalized by the total 
baryonic mass. Again, we compare two different injection 
efficiencies (£sn =0.1 and Csn = 0.3) with a reference 
case where no cosmic ray physics is included. In general, 
star formation is found to be most efficient at intermediate 
mass scales of ~ 1O 11 M0 in these simulations. However, 
the simulations with cosmic ray production show a clear 
reduction of their integrated star formation rate for halos 
with mass below a few times 10 11 /i~ 1 M Q , an effect that 
becomes progressively stronger towards lower mass scales. 
For the 10 9 /i _1 M Q halo, the suppression reaches more 
than an order of magnitude for (sn = 0-3. 

In Figure IT21 we take a closer look at the spatial dis- 
tribution of the cosmic ray pressure in the different cases, 
and the profiles of the stellar disks that form. To this end, 
we show the projected gas density distribution in an edge- 
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Fig. 12. Effect of cosmic ray feedback on star formation in simulations of isolated disk galaxy formation. Each row 
shows results for a different halo mass, for Mhaio = 10 9 , 10 10 , 10 11 and 10 12 h~ 1 M Q from top to bottom. We compare 
the projected gas density fields at time t = 2.0 Gyr of runs without cosmic ray feedback (left column) to that of runs 
with cosmic ray production by supernovae (middle column). The gas density field is colour-coded on a logarithmic 
scale. For the simulation with cosmic rays, we overplot contours that show the contribution of the projected cosmic ray 
energy density to the total projected energy density (i.e. thermal plus cosmic rays), with contour levels as indicated 
in the panels. Finally, the right column compares the azimuthally averaged stellar surface density profiles at time 
t = 2.0 Gyr for these runs. Results for simulations without cosmic ray physics are shown with a solid line, those for 
simulations with CR feedback with a dot-dashed line. 



22 



M. Jubelgas, V. Springel, T. Enfilin, C. Pfrommer: Cosmic rays in hydrodynamical simulations 



1.00 



5 





i i - 


no CR </ / / 


^ k ^-« - 






i i i i 





10 e 10 io 1Q n 10 H 10 i3 1Q 14 

Fig. 11. Efficiency of star formation as a function of halo 
mass in our isolated disk formation simulations. We show 
the ratio of the stellar mass formed to the total baryonic 
mass in each halo, at time t — 3.0 Gyr after the start of 
the simulations, and for two different efficiencies of cosmic 
ray production by supernovae. Comparison with the case 
without cosmic ray physics shows that star formation is 
strongly suppressed in small halos, by up to a factor ~ 
10 — 20, but large systems are essentially unaffected. 



on projection at time t = 2.0 Gyr, comparing the case 
without cosmic rays (left column) to the case with cos- 
mic rays (middle column), for a range of halo masses from 
1O 9 /i -1 M to lO 12 h^ 1 M . For the simulation with cos- 
mic rays, we overlay contours for the relative contribution 
of the projected cosmic ray energy to the total projected 
energy density. This illustrates, in particular, the spatial 
extent the cosmic ray pressure reaches relative to the star- 
forming region. Finally, the panels on the right compare 
surface density profiles of the stellar mass that has formed 
up to this time. 

Consistent with our earlier results, the stellar density 
profiles of the low mass halos show a significant suppres- 
sion when cosmic rays are included, while they are essen- 
tially unaffected in the high mass case. Interestingly, we 
also see that the gaseous disks in the low mass halos ap- 
pear to be "puffed up" by the additional pressure of the 
cosmic rays. It is remarkable that in the two low-mass 
systems there is substantial CR pressure found signifi- 
cantly above the star-forming regions, at densities much 
below the star formation threshold. This is despite the 
fact the acceleration of relativistic particles only occurs 
in star- forming regions of high density within the galactic 
disk in these simulations. Presumably, some of the CR- 
pressurized gas buoyantly rises from the star-forming disk 
into the halo, a process that is suppressed by the stronger 
gravitational field in the high mass systems. 

As a final analysis of our isolated disk simulations, we 
examine how well our simulation methodology for cosmic 
ray feedback converges when the numerical resolution is 
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Fig. 13. Resolution study of the star formation rate dur- 
ing the formation of a galactic disk in a halo of mass 
lO lo ft. _1 M , including production of CRs with an effi- 
ciency of £sn = 0.1. We compare results computed with 
10 5 , 4 x 10 5 , and 1.6 x 10 6 gas particles, respectively. 



varied. To this end we repeat one of the simulations with 
cosmic ray feedback (£sn = 0.1) of the 10 11 /i _1 Mq halo 
using a higher number of gas particles, namely 4 x 10 5 
and 1.6 x 10 6 , respectively. In Figure IT3l we compare the 
resulting star formation rates. While there are some small 
fluctuations when the resolution is varied, there is clearly 
no systematic trend with resolution, and the results ap- 
pear to be quite robust. In particular, the star formation 
rates for the simulations with 10 5 and 1.6 x 10 6 parti- 
cles are in very good agreement with each other despite 
a variation of the mass resolution by a factor of 16. Note 
also that the oscillations are reproduced by all three res- 
olutions, but they are not exactly in phase. Overall, this 
resolution test is very promising and suggests that the 
numerical model is well posed and can be applied to cos- 
mological simulations where the first generation of galax- 
ies is typically only poorly resolved. We can still expect 
meaningful results under these conditions. 

5.2. Cooling in isolated halos 

The comparison of the maximum cosmic ray cooling 
timescale with the thermal cooling time in the bottom 
panel of Figure [S] has shown that for a relatively wide tem- 
perature range, the lifetime of CRs is larger than the ther- 
mal cooling time. In a composite gas with a substantial 
cosmic ray pressure component, this could potentially sta- 
bilize the gas temporarily and reduce the rate at which gas 
cools and accumulates at the bottom of the potential well 
of a halo. Also, models have been conjectured where the 
temperature structure of the intracluster medium, with its 
characteristic observed minimum of one-third of the am- 
bient cluster temperature, could be explai ned by a str ong 
CR component in the intracluster medium ICerJ l)2005|) . 
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Fig. 14. Relative suppression of star formation in simula- 
tions of isolated halos when a fraction of 0.3 of the initial 
thermal pressure is replaced by a CR component of equal 
pressure. We show results as a function of halo virial mass 
for two different times after the simulations were started, 
for t = 1.0 Gyr (solid line), and for t = 3.0 Gyr (dot- 
dashed). For halos of low mass, the cosmic ray pressure 
contribution can delay the cooling in the halos. 

We here want to get an idea about the potential 
strength of this effect, and examine to this end a small 
set of toy simulations. To this end we consider a series 
of self-similar dark matter halos with a gas component 
that is initially in hydrostatic equilibrium, just as before 
in Section IH~T1 In fact, we use the same initial conditions 
as before, except that we replace a fraction of the initial 
thermal pressure with cosmic ray pressure, keeping the to- 
tal pressure constant. For definiteness, we choose a spec- 
tral cut-off of q = 1.685 and a spectral index a = 2.5 for 
the initial CR population. We then evolve the halos for- 
ward in time, including radiative cooling processes for the 
thermal gas as well as cosmic ray loss processes, but we 
disregard any sources of new cosmic ray populations. We 
are interested in the question whether the cooling flows 
that develop in these halos are modified by the presence 
of the cosmic rays. Note that some studies have predicted 
cosmic ray pressure contributions of up to ~ 50 per cent i n 
clusters of galaxies l|Miniati et allkOOlllRvu et alll2003|) . 
These fiducial test simulations can give a first indication 
of the size of the change of the cooling rates if these claims 
are indeed realistic. 

In Figure El we show the results of these simulations 
as a function of halo mass, in the form of the integrated 
star formation rate relative to an equivalent simulation 
without initial CR population. The cumulative star for- 
mation activity can here be taken as a proxy for the inte- 
grated strength of the cooling flow in the halo. We see 
that the total star formation for cluster halos of mass 
10 15 h~ 1 M & is essentially unchanged in the first 1 — 2 Gyr 
of evolution, while at late times, it is even slightly in- 



creased. For systems of significantly lower mass, the star 
formation rates are however reduced in the CR case, by 
up to ~ 40 per cent. This can be qualitatively understood 
based on a comparison of the thermal radiative cooling 
time with the CR dissipative cooling time. As the lower 
panel of Fig.JSJhas shown, the timescales are comparable at 
the virial temperature corresponding to the 10 15 /i _1 M© 
halo, but are quite different for lower temperatures, where 
the radiative cooling is significantly faster. In fact, a naive 
comparison of these timescales would perhaps suggest an 
even stronger suppression of the cooling efficiency in sys- 
tems of intermediate mass. In reality, the effect turns out 
to be much more moderate. This can be understood based 
on the softer equation of state of the CR component, com- 
bined with the fact that its cooling timescale typically de- 
clines faster than that of thermal gas when a composite 
gas is compressed. As a result, the ability of a CR pres- 
sure component to delay thermal collapse for a long time 
is quite limited, unless perhaps active sources for new pop- 
ulations of CR particles are present. 

6. Cosmological Simulations 

Previous simulation work on the effects of cosmic rays on 
structure formation has not accounted for the dynamical 
effects due to cosmic ray pressure, i.e. the effectiveness 
of cosmic ray production has only been estimated pas- 
sively. Interestingly though, these works suggested that 
the cosmic ray production may be quite ef ficient, with up 
to ~ 50% of the pressure being due to CRs llMiniati et all 
l200lURvu et alll2003tlRvu fc Kand.E)oil2004 . Here we 
present the first self-consistent cosmological simulations 
of CR production that also account for the dynamical ef- 
fects of cosmic rays. We first study the global efficiency 
of cosmic ray production at structure formation shocks. 
We then study the influence of cosmic ray feedback on 
star formation in cosmological simulations, and on possi- 
ble modifications of the Lyman-a forest. Finally, we study 
the modification of the thermodynamic properties of the 
intracluster medium in high-resolution simulations of the 
formation of a cluster of galaxies. 

6.1. Cosmic ray production in structure formation 
shocks 

In this subsection, we examine the efficiency of cosmic ray 
production at structure formation shock waves. To this 
end, we use simulations that include cosmic ray produc- 
tion injection at shocks and the cosmic ray loss processes 
(i.e. Coulomb cooling and hadronic losses), but we dis- 
regard radiative cooling and star formation. The cosmo- 
logical model we have simulated is a concordance ACDM 
model with parameters ACDM model, with parameters 
O = 0.3, ct 8 = 0.84, baryon density Q b = 0.04. We have 
picked a comoving box of side-length L = 100 h~~ 1 Mpc, 
and simulated each of our models at two resolutions, with 
2 x 128 3 and 2 x 256 3 particles, respectively. The results 
of the two resolutions arc in good agreement with each 
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other, so we restrict ourselves to reporting the results of 
the higher resolution runs with 2 x 256 3 particles in the 
following. 

We compare three simulations that differ in the treat- 
ment of the cosmic ray physics. In our 'full model', we 
account for shock injection self-consistently, i.e. we use 
the Mach number est imator of our companion paper 
^Pfrommer et all 1200^ to determine the energy content 
and the slope of the proton populations accelerated at each 
shock front. This simulation hence represents our best 
estimate for the overall efficiency of the CR production 
process due to structure formation shocks. We contrast 
this simulation with a model where the CR injection has 
been artificially maximized by adopting a fixed efficiency 
Cinj = 0.5 and a fixed injection slope a = 2.5 for all shocks. 
Note that this high efficiency is normally only reached as 
limiting case for high Mach number shocks, so that this 
model also allows us to assess the importance of the de- 
pendence of the shock injection efficiency on Mach num- 
ber. Finally, we compare these two models with a reference 
simulation where no cosmic ray physics was included. This 
reference simulation is hence a standard non-radiative cal- 
culation where only shock-heating is included and the gas 
behaves adiabatically otherwise. 

In Figure 1151 we compare the time evolution of the 
mean mass-weighted temperature of the full cosmic ray 
model to that in the ordinary non-radiative simulation. 
We also include a measurement of the mean energy in 
cosmic rays, converted to a fiducial temperature using the 
same factor that converts thermal energy per unit mass to 
temperature. Interestingly, at high redshift the cosmic ray 
energy content evolves nearly in parallel to the thermal 
energy, and both are roughly half what is obtained in the 
simulation without cosmic rays. Apparently, the thermal- 
ization of gas is dominated by strong shocks which reach 
the asymptotic injection efficiency of 50 percent. At late 
times, however, the CR energy does not grow as quickly 
as the thermal energy content any more, and the thermal 
energy in the CR simulation becomes closer to the thermal 
energy in the ordinary simulation. 

These trends become more explicit when the energy 
content in CRs and in the thermal reservoir of the full 
CR simulation is divided by the thermal energy content 
of the reference simulation, as shown in the bottom panel 
of Figure ^] Around rcdshifts z ~ 6 — 10, the CR en- 
ergy content nearly reaches the same value as the thermal 
energy in the full CR-model, and their sum is essentially 
equal to the thermal energy in the simulation without cos- 
mic rays. Over time, the relative importance of the cosmic 
rays declines, however, and the thermal energy in the full 
CR model slowly climbs back to the value obtained in the 
non-radiative reference simulation. At the same time, the 
sum of cosmic ray and thermal energy obtained in the full 
model becomes a few percent higher at the end than that 
in the simulation without cosmic rays, despite the fact 
that the CR simulation loses some energy to radiation via 
the hadronic decay channels. Apparently, the simulation 
with cosmic rays extracts slightly more energy out of the 
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Fig. 15. Time evolution of the mean thermal energy and 
the cosmic ray energy content of the gas in non-radiative 
cosmological simulations. In the top panel, the solid thick 
line shows the mass- weighted temperature for a simulation 
where the efficiency of cosmic ray production at struc- 
ture formation shocks is treated self-consistently based on 
our on-the-fly Mach number estimator. The dashed line is 
the corresponding mean cosmic ray energy, which we here 
converted to a fiducial mean temperature by applying the 
same factor that converts thermal energy per unit mass 
to temperature. For comparison, the thin solid line shows 
the evolution of the mean mass-weighted temperature in 
an ordinary non-radiative simulation without cosmic ray 
physics. In the bottom panel, we show the ratio of the 
mean thermal energy in the cosmic ray case relative to 
the energy in the simulation without cosmic rays (solid 
line), while the dashed line is the corresponding ratio for 
the cosmic ray component. Finally, the dotted line gives 
the total energy in the cosmic ray simulation relative to 
the ordinary simulation without cosmic rays. 



gravitational potential wells of the dark matter. An expla- 
nation for this behaviour could derive from the fact that 
more energy needs to be invested into CRs to reach the 
same pressure compared with ordinary thermal gas. This 
should allow the gas in CR simulations to fall deeper into 
gravitational potential wells before it is stopped by shocks 



M. Jubelgas, V. Springel, T. Enfilin, C. Pfrommer: Cosmic rays in hydrodynamical simulations 



25 




Fig. 16. Mean relative contribution of the cosmic ray pressure to the total pressure, as a function of gas overdensity 
in non-radiative cosmological simulations. We show measurements at epochs z = Q, 1, 3, and 6. The panel on the 
left gives our result for a simulation where the injection efficiency and slope of the injected cosmic ray spectrum are 
determined based on our on-the-ffy Mach number estimation scheme. For comparison, the panel on the right shows 
the result for a simulation with a fixed injection efficiency Clin = 0.5 and a soft spectral injection index of a- m ] = 2.9. 
Clearly, the relative contribution of cosmic ray pressure becomes progressively smaller towards high densities. It is 
interesting that the trends with redshift are reversed in the two cases. 



and pressure forces, such that more gravitational energy 
is liberated overall. 

It is also interesting to ask how the relative impor- 
tance of cosmic rays depends on gas density. This is ad- 
dressed by Figure 1161 where we show the relative con- 
tribution of cosmic rays to the total gas pressure, as a 
function of baryonic overdensity, separately for different 
redshifts. We give results both for the simulation with 
self-consistent shock injection (left panel), and for the one 
where we imposed a constant injection efficiency (right 
panel). In general, the importance of cosmic rays is largest 
for densities around the mean cosmic density, and de- 
clines towards higher densities. This is consistent with 
the expectation that the strongest shocks occur at low 
to moderate overde nsities in the accretion regions around 
halos and filaments faang et al. , ll996HQuilis et all Hi 
iMiniati l2002t iRvu et all l2003t IPfrommer et all 12006 



and also with the growing importance of cosmic ray 
loss processes at high densities. Another interesting trend 
found in our self-consistent simulation (the 'full model) is 
that cosmic rays are particularly important at high red- 
shift, with a gradual decline towards lower redshift, sug- 
gesting that the mean Mach number of shocks becomes 
lower at later times, as is indeed confir med by studies of 
the cosmic Mach num ber distribution I Rvu et all [2003: 
IPfrommer et all l200fit) . Note that this trend is reversed 
in our fiducial simulation with a fixed shock injection effi- 
ciency, where at all overdensities the relative importance 
of CRs grows with cosmic time. This emphasizes that a 
correct accounting of the shock strengths is highly im- 
portant even at a qualitative level to correctly model the 
evolution of the cosmic ray pressure distribution. We note 



that an implicit assumption we made in the above analysis 
is that weak magnetic fields are ubiquitous in the universe, 
even at low density. Whether they really exist and where 
they ultimately come from is an open question however. 

In Figure El we show the projected gas density field 
in a slice through the simulation box at z — 0. To high- 
light the relative importance of cosmic rays, the panel 
on the right shows the ratio of the projected cosmic ray 
energy density to the projected thermal energy density. 
The relative importance of CRs is clearly highest in the 
volume-filling gas at low density. In the accretion regions 
around halos and filaments, the CR contribution is still 
comparatively large, but the high-density regions inside 
massive halos are avoided, in agreement with the results of 
Figure^] This raises the interesting question whether cos- 
mic rays may perhaps modify the state of the intergalactic 
medium to the extent that the properties of Lyman-a ab- 
sorption systems are modified. The latter arise primarily 
in gas that is largely unshocked, so that the effects might 
be weak even though the CR pressure contributions are 
predicted to be on average large exactly at overdensities 
of a few. We will come back to this question in Sect ion RT3l 

Another interesting question is whether the bulk prop- 
erties of halos are modified by the CR production at large- 
scale structure shocks. We are for example interested in 
the question whether the concentration of gas in halos 
is changed, which could manifest itself in a modification 
of the mean gas mass inside dark matter halos. To ex- 
amine this question we determine halo catalogues for our 
simulations by means of the FOF algorithm with a stan- 
dard linking- length of 0.2, and measure the virial radii and 
masses by means of the spherical overdensity algorithm. In 
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Fig. 17. Projected gas density field (left panel) in a slice of thickness 20 /i _1 Mpc through a non-radiative cosmological 
simulation at z — 0. The simulation includes cosmic ray production at structure formation shocks using self-consistent 
efficiencies based on an on-the-fly Mach number estimation scheme. The panel on the right shows the ratio of the 
projected cosmic ray energy density to the projected thermal energy density. We clearly see that the local contribution 
of cosmic rays is largest in voids. It is also still large in the accretion regions around halos and filaments, but is lower 
deep inside virialized objects. 



Figure^] we show the mean baryonic mass fraction in ha- 
los as a function of halo mass for the simulation with self- 
consistent CR injection and for the run without cosmic 
ray physics. In both cases, the baryonic fraction within the 
virial radius lies slightly below the universal baryon frac- 
tion, reaching ~ 0.91 — 0.94 of it, and for poorly resolved 
halos it drops a bit further. Such a depression of the uni- 
versal baryon fractio n is generally foun d in non-radiative 
SPH simulations fe.g lFrenk et allll999(l . However, a com- 
parison of the two simulations shows that the halos in the 
run with CRs have systematically increased their baryonic 
fraction, albeit by only about 1-2 per cent of the univer- 
sal baryon fraction. This is consistent with expectations 
based on the higher compressibility of a composite gas 
with thermal and CR components. 

Using the group catalogues, we can also measure the 
mean cosmic ray energy content inside the virial radius 
of halos. In Figure we show the ratio of cosmic ray to 
thermal energy as a function of halo mass. For the simula- 
tion with a fiducial shock injection efficiency of £j n j = 0.5 
at all shocks, the ratio we obtain is ~ 0.2, independent 
of halo mass. Loss processes in the CR population and 
the shallower adiabatic index of the CR component make 
this value much smaller than expected in this case for the 
post-shock region of a single shock where one would ex- 
pect ~ 0.5. In the simulation with self-consistent injection 
of CRs, we find an interesting mass dependence where 
the ratio of CR-to-thermal energy drops from about 0.2 
for 10 12 /i-'Mn halos to ~ 0.05 for 10 15 h^Mp, halos. 



Apparently for building up the thermal energy of clus- 
ters of galaxies, weaker shocks and adiabatic compression 
are comparatively more important than for galaxy-sized 
halos. We note that the value of ~ 5 — 10% we predict 
here for the CR energy content due to structure forma- 
tion shocks in clus ters of galaxies is qui te a bit lower than 
previous estimates iMiniati et al.l (|200lf ) . However, it is in 
good agreement wit h CR constraints from gamma ra y and 
radio observations IjPfrommer fc En filin. 2003, 2004). 

6.2. Dwarf galaxy formation 

We now turn to studying the effects of cosmic ray feedback 
on galaxy formation in cosmological simulations. We have 
already found that small galaxies should be affected most. 
We hence expect small dwarf galaxies to be most suscepti- 
ble to sizable effects of CR feedback from star formation. 
To reach a reasonably good mass resolution, we simulate 
periodic boxes of side-length 10ft. _1 Mpc, using 2 x 256 3 
particles. This gives a mass resolution of 6.62 x 10 5 A^'Mq 
and 4.30 x 10 6 /i _1 Mq in the gas and dark matter, respec- 
tively. We limit ourselves to evolving the simulations to 
a redshift of z = 3, because at lower redshift the funda- 
mental mode of the small simulation volume would start 
to evolve non-linearly, at which point the simulation as 
a whole could not be taken as representative for the uni- 
verse any more. We are hence restricted to studying the 
high-rcdshift regime, but we expect that our results are 
indicative for the trends that would be seen in the dwarf 




Fig. 18. Mean baryon fraction within the virial radius as a 
function of halo mass, normalized by the universal baryon 
fraction. We compare results for two non-radiative sim- 
ulations, one with cosmic ray production by shocks, the 
other without cosmic ray physics. The bars show the la 
scatter among the halos in each bin. When cosmic rays are 
included, the compressibility of the gas in halos becomes 
larger, leading to a slight increase of the enclosed baryon 
fraction. 

galaxy population at lower redshifts as well, provided suf- 
ficiently well resolved simulations are available. 

We have simulated the same initial conditions three 
times, varying the cosmic ray physics included. The first 
simulation is a reference run where we only included radia- 
tive cooling and star formation but no cosmic ray physics. 
The second simulation is a model where we also consid- 
ered cosmic ray production by the supernovae associated 
with star formation, using an efficiency of Csn = 0.35. 
Finally, our third simulation is a model where we in addi- 
tion included cosmic ray production by structure forma- 
tion shocks, using the self-consistent efficiencies derived 
from our on-the-fly Mach number estimation scheme. The 
latter simulation hence represents our best estimate for 
the total effect of cosmic rays on dwarf galaxy formation. 

In Figure 1201 we compare the cosmic star formation 
histories of the three simulations up to a redshift of 
z ~ 2.9. The incorporation of cosmic ray production by 
supernovae leads to a significant reduction of the high red- 
shift star formation activity, but the shape of the star for- 
mation history, in particular its exponential rise, is not 
changed significantly. At these high redshifts, the star for- 
mation rate is dominated by small halos which are affected 
strongly by CR feedback, so this result is not unexpected 
given our previous findings. If CR production by structure 
formation shocks is included as well, the star formation is 
reduced further, although by only a small factor. This in- 
dicates that the cosmic ray pressure component injected 
into forming halos indeed tends to slightly slow down the 



Fig. 19. Ratio of energy in cosmic rays to thermal energy 
within the virialized region of halos, shown as a func- 
tion of halo mass. We compare results for two different 
non-radiative simulations, one treating the production of 
cosmic ray at shocks fronts using a self-consistent Mach 
number estimator, the other invoking a constant injec- 
tion efficiency. The bars give the la scatter among the 
halos in each bin. Interestingly, the self-consistent injec- 
tion scheme predicts a lower CR energy content in more 
massive systems. In contrast, a constant shock injection 
efficiency produced no significant trend of the CR energy 
content with halo mass. 

cooling rates, consistent with the results we found for iso- 
lated halos. Towards redshift z ~ 3, the differences in the 
star formation rates become noticeably smaller however, 
suggesting that the low redshift star formation histories 
will differ at most by a small amount. Since the bulk of 
the star formation shifts to ever larger mass scales at low 
redshift ijSpringel fc Hernauisd 12003b]) . this can be easily 
understood in terms of the smaller influence of CR feed- 
back on large halos. 

In order to make the effects of CR feedback on small 
halos more explicit, we determine halo catalogues in the 
simulations using a group finder. We are especially inter- 
ested in the question how the efficiency of star formation 
is changed by the inclusion of cosmic rays as a function of 
halo mass. In Figure l2"Tl we show the total-to-stellar mass 
ratios of these groups as a function of halo mass, both for 
the simulation with CR production by supernovae, and 
for the simulation without cosmic ray feedback. The sim- 
ulation where also CR production by shocks is included 
is quite similar on this plot to the simulation that only 
accounts for CR from supernovae, and is therefore not 
shown. The symbols show the mean total-to-stellar mass 
ratio in small logarithmic mass bins, while the bars in- 
dicate the scatter by marking the central 68% percentile 
of the distribution of individual halos. The results show 
clearly that the star formation is significantly reduced by 
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Fig. 20. Evolution of the cosmic star formation rate den- 
sity in simulations of galaxy formation at high redshift. 
We compare results for three simulations that include dif- 
ferent physics, a reference simulation without cosmic ray 
physics, a simulation with CR production by supernovae, 
and a third simulation which in addition accounts for CR 
acceleration at structure formation shocks with an effi- 
ciency that depends on the local Mach number. 

CRs for low-mass halos, by factors of up to « 10 for host 
halo masses of ~ 10 10 Mq/i _1 and below. On the other 
hand, the amount of stars produced in massive halos is 
hardly changed. It is particularly interesting that the ef- 
fect of CRs manifests itself in a gradual rise of the total- 
to-stellar mass ratio towards lower masses. This can be 
interpreted as a prediction for a steeply rising 'mass-to- 
light' ratio towards small halo masses, which is exactly 
what appears to be needed to explain the observed lu- 
minosity function of galaxies in the ACDM concordance 
model. The problem is here that the halo mass function 
increases steeply towards low mass scales. If the mass-to- 
light ratio is approximately constant for low masses, this 
leads to a steeply rising faint end of the galaxy luminosity 
function, in conflict with observations. However, a steeply 
rising mean mass-to-light ratio towards low mass halos 
could resolve this problem and provide a suitable 'trans- 
lation' between the halo mass function and the galaxy 
luminosity function. 

We note that conditional luminosity function analysis 
of the 2 D egree Field Galaxy Redshift Survey (2dFGRS) 
has shown l|van den Bosch et alll2003|) that there appears 
to be a minimum in the observed mass-to-light ratio of 
galaxies around a halo mass of sw 3x 10 11 h M®. This 
feature is reproduced surprisingly well in our simulations, 
although even with CR feedback included, the rise of the 
stellar mass to light ratio towards low masses appears to 
be not as sharp as required based on their analysis. 

However, one needs to caution that the results of 
Fig. |^ cannot be naively translated into changes of the 
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Fig. 21. Comparison of the averaged total mass-to-light- 
ratio within the virial radius of halos formed in two high- 
resolution cosmological simulations up to z = 3. Both sim- 
ulations follow radiative cooling and star formation, but 
one also includes CR-feedback in the form of cosmic pro- 
duction by supernovae, with an efficiency of £sn = 0.35 
and an injection slope of q;sn = 2.4. The bars indicate the 
scatter among halos in the logarithmic mass bins (68% 
of the objects lie within the range marked by the bars). 
Clearly, for halo masses below 10 11 h~ 1 M Q , CR feedback 
progressively reduces the overall star formation efficiency 
in the halos. 



faint-end slope of the luminosity function, as seen when 
we directly compare the K-band luminosity functions at 
z — 3. To determine those, we identify individual groups 
of stars as galaxies using a mod ification of the SUBFIND 
algorithm ijsnringel et a.lll200lft for detecting bound sub- 
structures in halos. For each of the galaxies, we esti- 
mate magnitudes in stand ard observational band based on 
iBruzual fc Charlotl l|2003f) population synthesis models. In 
Figure [221 we compare the resulting restframe K-band lu- 
minosity functions at z — 3 for the simulations with CR 
feedback by supernovae and the simulation without any 
cosmic ray physics. We see that both luminosity functions 
are well fit by Schechter functions, with faint-end slopes 
of a = 1.15 and a = 1.10, respectively, for the cases with- 
out and with CR feedback. We hence find that CRs only 
mildly reduce the faint-end slope despite their differen- 
tial reduction of the star formation efficiency towards low 
mass scales. The result needs to be taken with a grain of 
salt though, as the faint-end slope could still be influenced 
by resolution effects in these simulations. A final assess- 
ment of the importance of CR feedback in shaping the 
faint-end of the galaxy luminosity function needs there- 
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Fig. 22. The K-band galaxy luminosity functions at z = 3 
in two high-resolution cosmological simulations. One of 
the simulations follows ordinary radiative cooling and star 
formation only (blue), the other additionally includes cos- 
mic ray production by supernovae (red). The latter re- 
duces the faint-end slope of the Schechter function fit 
(solid lines) to the data measured from the simulations 
(histograms). It is reduced from —1.15 to —1.10 in this 
case. 



fore await future simulations with substantially increased 
resolution. 



6.3. Cosmic ray effects on the intergalactic medium 

As the Mach number distribution is dominated by strong 
shocks at high redshift, we expect that cosmic ray pro- 
duction is particularly efficient at early epochs and at the 
comparatively low densities where the strongest shocks oc- 
cur, provided sufficient magnetization of the IGM existed 
to allow CR acceleration to operate. Also, the thermaliza- 
tion time scales of cosmic rays are quite long at low den- 
sities. Figure UHl has shown that the mean energy content 
of cosmic rays can reach a sizable fraction of the thermal 
energy content at around redshift z ~ 3, suggesting a po- 
tentially important influence on the intergalactic medium 
at this epoch. Note however that in computing the results 
of Figure 1161 we had neglected cosmic reionization, which 
will boost the thermal energy relative to the cosmic ray 
content. Also, large parts of the IGM at z = 3, particularly 
those responsible for the absorption seen in the Lyman-a 
forest, consist largely of unshocked material. Whether the 
Lyman-a forest might show any trace of the influence of 
cosmic rays is therefore an interesting and open question. 

To investigate this question further, we have com- 
puted Ly-a absorption spectra for the cosmological sim- 



0.1000 



0.0100 



0.0010 



0.0001 




z = 3 



with cosmic rays 
without cosmic rays 
Mc Donald et al. (2000) 



0.01 



0.10 

k [ km 1 s ] 



1.00 




0.01 



0.10 



1.00 



k [ km 1 s ] 



Fig. 23. Ly-a flux power spectrum (top) at z — 3 in sim- 
ulations with and without cosmic ray production in struc- 
ture formation shocks. The results lie essentially on top 
of each other, and only by plotting their ratio (bottom 
panel), it is revealed that there are small differences. In 
the simulation with cosmic rays, the power is suppressed 
by up to ~ 15% on scales 0.1km _1 s < k < 0.7km _1 s, 
while there is an excess on still smaller scales. However, 
on large scales k < 0.1 km -1 s, which are the most rel- 
evant for determinations of the matter power spectrum 
from the Ly-a forest, the power spectrum is not changed 
by including CR physics. For c omparison, we have also in- 
cluded observational data from iMcDonald et all (l200Ct) m 
the top panel (the open symbols are corrected by remov- 
ing metal lines) . A slightly warmer IGM in the simulations 
could account for the steeper thermal cut-off observed in 
the data. 



ulations with 10/i _1 Mpc boxes analysed in the previous 
section. The two simulations we have picked both include 
radiative cooling, star formation, and heating by a spa- 
tiall y uniform UV backgrou nd bases on a slightly modi- 
fied [HaaKl^^^^adau| l|l996|) model, with reionization at 
redshift z = 6. While one of the simulations did not ac- 
count for any cosmic ray physics, the other included cos- 
mic ray production by large-scale structure shocks and 
supernovae, as well as dissipative loss processes in the CR 
population. 

For both simulations, we computed Lyman-a absorp- 
tion spectra for 2048 lines of sights, along random di- 
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rections parallel to the principal axes of the simulation 
boxes. By slightly adjusting the UV intensity, we have 
renormalized the spectra to the same mean transmission 
of (t) = 0.68. A direct comparison of the spectra along 
the same lines-of-sight through the two simulations shows 
essentially perfect agreement, with very small residuals. 
This already indicates that any systematic difference be- 
tween the simulations must be quite subtle, if present. To 
objectively quantify this, we have computed the average 
1-d flux power spectra for the two cases and compare them 
in Figure E3I The top panel compares the two flux spec- 
tra directly with each o ther, and to observational data of 
iMcDonald et alJ l)2000|) . The results for the two simula- 
tions lie essentially on top of each other in this represen- 
tation. The agreement with observational data is good, 
apart from a small excess of power on small scales, which 
can however be understood as a consequence of the too 
cool temperature of the IGM in our simulations compared 
with observations. 

More interesting is perhaps an examination of the 
ratio of the flux power spectrum with cosmic rays to 
that without cosmic rays, as shown in the bottom panel 
of Figure [531 While for large-scale modes with k < 
0.1km _1 s, no noticeable differences are seen, there is 
a 5-15% reduction of power in the wave-length range 
0.1 km _1 s < k < 0.7km _1 s, and at still smaller scales, the 
difference changes sign and turns into a growing excess of 
power in the CR simulation. These effects of CRs on the 
Ly-a therefore lie in a regime that is normally not used 
to constrain the matter power spectrum with Lyman-a 
forest data, at least in conservative t reatments that focus 
on k < 0.03 km^s (IViel et all Eool. In general we hence 
find that the effects on the Lyman-a forest are very small 
and subtle; the forest survives CR injection by large-scale 
structure shocks essentially unaltered, even though they 
contribute a sizable fraction to the mean energy content of 
the gas due to shock dissipation at densities at and around 
the mean density of the universe. Note that our simula- 
tions did not allow for a possible diffusion of CRs, but it 
seems unlikely that including this effect could change this 
conclusion. 

6.4. Formation of clusters of galaxies 

In this section, we study in more detail the influence of cos- 
mic rays on individual halos formed in cosmological simu- 
lations. We focus on high-resolution 'zoom' simulations of 
the formation of a massive cluster of galaxies. Such 'zoom' 
simulations are resimulations of an object identified in 
a cosmological structure formation simulation with large 
box-size. Once the object of interest has been selected, its 
particles' are traced back through time to their origin in 
the unperturbed initial conditions. The Lagrangian region 
of the cluster thus identified is then populated with many 
more particles of lower mass, thereby increasing the local 
resolution, while in regions further away, the resolution is 
progressively degraded by using ever more massive parti- 
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Fig. 24. Spherically averaged radial profiles of thermody- 
namic gas properties in three re-simulations of the same 
cluster of galaxies. All three simulations were not follow- 
ing radiative cooling and star formation, and the reference 
simulation shown with a solid line does not include any 
CR physics. However, the simulation shown with a dashed 
line accounted for CR production at structure formation 
shocks with a fixed efficiency (Cinj = 0.5, ai n j = 2.9) while 
for the simulation shown with dot-dashed lines, the shock 
injection efficiency was determined self-consistently based 
on our Mach number estimation scheme. The panel on top 
compares the thermal pressure in the three simulations. 
For the two simulations with cosmic rays, we additionally 
plot the CR-pressure, marked with symbols. The panel in 
the middle compares the gas temperature, while the panel 
on the bottom shows the radial run of the gas density. The 
vertical dotted line marks the virial radius of the cluster. 



M. Jubelgas, V. Springel, T. Enfilin, C. Pfrommer: Cosmic rays in hydrodynamical simulations 



31 



cles. In this way, the computational effort can be concen- 
trated in the object of interest, while at the same time its 
cosmological environment is still accounted for accurately 
during its formation. 

We have computed 6 resimulations of the same clus- 
ter of galaxies, using different models for the physics of 
radiative cooling, star formation, and cosmic rays. The 
cluster has been selected from a set of zoome d cosmologi- 
cal init ial conditions originally constructed bv lDolag et al.l 
I 2004a) and has a virial mass of wlO 14 ft _1 M0 at red- 
shift z = 0. The gas particle mass is 1.6x 10 s /i~ 1 M Q in 
the high resolution region, implying that the cluster is re- 
solved with roughly 300000 gas and 300000 dark matter 
particles within the virial radius. The gravitational soften- 
ing length for the simulations was kept fixed in comoving 
units at redshifts z > 5, and then held constant in physical 
units at a value of 5 /i~ 1 kpc at lower redshifts. 

Our 6 simulations fall into two groups of 3 each. In 
the first group, we have not included radiative cooling 
processes and star formation. Here we use a non-radiative 
('adiabatic') simulation as a reference run, and compare it 
with two simulations that include cosmic ray production 
at large-scale structure formation shocks, one of them us- 
ing the self-consistent Mach-number dependent injection 
efficiency, and the other a fixed efficiency of Cinj = 0.5 with 
ctinj = 2.5 for all shocks. This set hence parallels the types 
of simulations analyzed in section liTTl In the second set of 
3 simulations, we include radiative cooling and star forma- 
tion. Again, we consider one reference simulation without 
any cosmic ray physics, and compare it with two simula- 
tions where cosmic rays are included. These two consist 
of one run where cosmic rays are only injected by super- 
novae associated with star formation (using £sn = 0.35, 
and «sn = 2.4), while the other also accounts for cosmic 
rays produced at shock waves. The second set of simula- 
tions hence corresponds to the types of simulations ana- 
lyzed in section l?T2l In all simulations with cosmic rays, we 
have assumed and index a = 2.5 and included Coulomb 
cooling and hadronic losses for the CR populations. 

In Figure 1241 we compare spherically averaged radial 
profiles of pressure, temperature, and gas density for the 
three non-radiative simulations. For the pressure, we show 
the thermal as well as the cosmic ray pressure for the 
two runs with cosmic ray physics. Interestingly, the cosmic 
ray pressure component is substantially below the thermal 
one, even in the fiducial case of a constant shock injection 
efficiency of Cinj = 0.5 for all shocks. However, in this case 
the thermal pressure is substantially elevated compared 
to the run without cosmic rays. This goes along with an 
increase of the gas density in the inner parts, and a re- 
duction of the thermal temperature throughout the clus- 
ter volume. This is the expected behaviour based on the 
higher compressibility of the gas in this case. 

However, the cosmic ray pressure in the simulation 
with a self-consistent injection efficiency is substantially 
smaller, and even at the virial radius is at most ~ 10 per- 
cent of the thermal pressure, while for much of the inner 
parts, r < 100ft. _1 kpc, the cosmic ray pressure contri- 
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Fig. 25. Spherically averaged radial profiles of thermody- 
namic gas properties in three re-simulations of the same 
cluster of galaxies. All three simulations included radiative 
cooling of the gas, star formation and supernova feedback. 
The solid lines show the results of a reference simulation 
which did not include any cosmic ray physics. The other 
two simulations included CR production by supernovae, 
and the one shown with dot-dashed lines in addition ac- 
counted for CR injection at structure formation shocks, 
using self-consistent efficiencies based on our Mach num- 
ber estimation scheme. The panel on top compares the 
thermal pressure in the three simulations. For the two 
simulations with cosmic rays, we additionally plot the CR- 
pressure marked with symbols. The panel in the middle 
compares the gas temperature for the three cases, and the 
panel on the bottom shows the radial run of the gas den- 
sity. 
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Fig. 26. Cumulative radial stellar mass profile in three 
re-simulations of the same cluster of galaxies. The simula- 
tions are the same ones also shown in FigureEHl The solid 
line gives the result for a reference simulation without CR 
physics, the dashed line includes CR production by super- 
novae, and the dot-dashed line additionally accounts for 
CR injection at structure formation shocks. The vertical 
dotted line marks the virial radius of the cluster. 



bution drops to the percent level and below. Obviously, 
cosmic rays are not produced efficiently enough at shocks 
to fill much of the ICM with a dynamically significant cos- 
mic ray pressure component, consistent with the results of 
Figure El Consequently we find that in this case the pro- 
files of gas density, temperature and thermal pressure are 
very similar to the corresponding results for the simulation 
without cosmic ray physics. 

In Figure we show the equivalent results for the 
radial profiles for the 3 simulations that include radiative 
cooling and star formation. Compared to the non-radiative 
calculations, the ICM shows a markedly different struc- 
ture. Due to the presence of a strong cooling flow, the 
temperature profile rises towards the centre, until it even- 
tually drops sharply at around 20 kpc due to the onset of 
efficient cooling. The gas density has become significantly 
lower in the bulk of the cluster volume due to the large 
amount of gas that has cooled out, and correspondingly, 
the total pressure has fallen in much of the cluster vol- 
ume. But there are also interesting differences in the sim- 
ulations with and without cosmic rays. Recall that both 
simulations included cosmic ray production by supernovae 
feedback, while only one of them accounted for cosmic rays 
by structure formation shocks. The CR pressure contribu- 
tion in both simulations is quite similar through most of 
the cluster, at the level of a few percent of the thermal 
pressure. Also note that in the very inner parts, where 
the gas drops out through cooling, the cosmic ray pres- 
sure rises sharply, even reaching and exceeding the ther- 
mal pressure. In this small region, the thermal pressure is 
dissipated more rapidly than the cosmic ray pressure. 



Finally, in Figure[2n|we compare the cumulative stellar 
profile of the cluster in the three simulations with radiative 
cooling and star formation. While the total mass of stars 
formed within the virial radius has become smaller by the 
inclusion of cosmic ray feedback, the stellar mass in the 
central cluster galaxy has actually increased. The cluster 
cooling flow has therefore slightly increased in strength, 
consistent with the results we obtained for isolated halos 
of this mass range. On the other hand, the luminosity of 
the smaller galaxies orbiting within the cluster has become 
smaller, in line with our finding that small galaxies expe- 
rience a reduction of their star formation activity. It is 
clear however that our results do not suggest that cosmic 
rays could provide a solution to the cooling flow problem 
in clusters of galaxies, at least not with the CR sources 
we have considered here. This conclusion could potentially 
change in interesting ways when CR pro duction by AGN 
in clusters of galaxies is inc luded as well I Churazov et all 
l2nmUF,nfllin fc Vogd . Eftf)^ . 

6.5. The influence of cosmic ray diffusion 

In all of our previous results, we have ignored the effects of 
cosmic ray diffusion, largely because of the uncertainty in- 
volved in constraining an appropriate diffusivity. However, 
diffusion could potentially be important in several envi- 
ronments, depending of course on the details of the mag- 
netic field structure and the strength of the resulting dif- 
fusivity. While our present formalism implemented in the 
simulation code is capable of dealing with isotropic dif- 
fusion, in reality the diffusion is likely to be anisotropic, 
governed by the local magnetic field configuration. In prin- 
ciple, cosmological structure formation calculations with 
SPH are capable of following magneto-hydrodynamic s 



llDolag et allll99flU200fiUPrice fc Monaghanll2004ll2005l) . 
although this is presently still fraught with numerical 
and physical difficulties. We therefore postpone a detailed 
analysis of the influence of cosmic ray diffusion to future 
work. Here, we investigate instead a simple example, that 
gives a first illustration of the effects that can be expected. 

To this end, we repeat our simulations of isolated disk 
galaxy formation with CR injection by supernovae, but 
this time with diffusion included. We use a parameter- 
ized diffusivity as described in section [21 setting the val- 
ues of the density and temperature scaling exponents to 
tit = 1/6 and n p = —1/2, respectively, with a base- 
line diffusivity of ~ 10kpc 2 Gyr^ 1 at the threshold for 
star formation, i.e. our diffusivity model is given by equa- 
tion Q54)). The simulations we repeat are the ones consid- 
ered in Section lETI with an injection efficiency of £sn = 0.3 
for the production of CRs by supernovae. 

In Figure E3 we compare the resulting star formation 
rates for halos of mass 10 9 ft _1 M0 and 10 10 /i _1 M Q as 
a function of time with the corresponding results without 
diffusion. Interestingly, the oscillations due to the unstable 
dynamics of a cosmic ray dominated ISM are substantially 
suppressed when diffusion is included. This effect is quite 
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Fig. 27. Effects of cosmic ray diffusion on the star formation and the pressure distribution in isolated halos of mass 
10 9 h~ 1 M.Q and 10 10 h~ 1 M & . The panels on top compare the star formation rate when CR diffusion is included (thick 
blue line) to the case where it is neglected (thin green line). The dotted lines show the result when CR-production 
by supernovae is not included. In the bottom panels, we show projected gas density fields through the halos at time 
t = 2.0 Gyr, with contours overlaid that give the fractional contribution of the projected CR energy to the total 
projected energy. These panels correspond directly to equivalent maps shown in Figure 1121 for the case without CR 
diffusion. 



strong in the results for the 10 9 /i _1 M Q halo, where we 
now observe a nearly constant, quiescent star formation 
rate. For the 10 10 /("'Mq halo, the oscillations are only 
partially washed out and happen less frequently, but if 
they occur, they are stronger. Here the star formation rate 
of the galaxy develops a 'bursty' character. Interestingly, 
diffusion actually reduces the integrated star formation 
still further; it drops by about 30% for the 10 9 /i _1 M 
halo, and by 21% for the 10 10 /! _1 Mq halo compared to the 
case without diffusion. Apparently, the cosmic rays that 
escape from the star-forming ISM into the less-dense gas 



in the halo are able to supply a partial pressure support 
that effectively reduces the rate at which gas cools. 

The more extended and smoother spatial distribution 
of cosmic rays due to diffusion can also be appreciated in 
the bottom panels of Figure l27l where we show projections 
of the gas density field with contours for the cosmic ray 
to thermal energy content overlaid. These panels directly 
correspond to the ones shown in Figure ^] for the case 
without diffusion. However, for halos of mass 10 11 /i _1 M© 
and more, diffusion with the strength considered here has 
no significant impact on the dynamics. The progressively 
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larger size of more massive systems makes it ever more 
important for diffusion to efficiently transport CR energy 
across the system. 

7. Conclusions 

In this paper, we have presented the details of the first 
practical implementation of a simulation code capable of 
carrying out high-resolution simulations of cosmological 
structure formation with a self-consistent treatment of 
cosmic ray physics. In particular, our method takes dy- 
namical effects of pressure forces due to cosmic rays into 
account and therefore allows us to carry out studies of CR 
feedback in the context of galaxy formation. The underly- 
ing formalism for the tre atment of cosmic ra ys is discussed 
m a companion paper (|Enfflin et all l2006(l and forms a 
compromise between the complexity of cosmic ray physics 
and the requirements of computational efficiency. In par- 
ticular, we use a simplified spectral representation in terms 
of power laws for the momentum distribution with a low 
momentum cut-off. This allows for a rather significant sim- 
plificationat the p rize of a moderate loss of accuracy. As 
Enfi lin et alJ l)2006|) have shown, the cosmic ray pressure 
is expected to be accurate to about 10 per cent in our 
model under steady state conditions. We argue that this 
is sufficiently accurate for our purposes given the other 
uncertainties and approximations involved. 

Our formalism also makes used of an on-the-fly shock 
detection scheme for SPH develop ed in a second compan- 
ion paper jPfrommer et alll2006|) . This method allows us 
to estimate Mach numbers in shocks captured during SPH 
simulations, such that we can use an appropriate efficiency 
for the CR injection at large-scale structure shock waves. 

We have given an initial analysis of the principal effects 
of two sources of cosmic rays, namely CRs produced by 
supernova explosions, and by diffusive shock acceleration 
during structure formation. The loss processes we con- 
sidered were Coloumb cooling and hadronic losses, which 
should be the most important ones. If desired, the mod- 
elling of these CR sink terms can be refined in the future 
within our methodology, and additional sources like cos- 
mic rays from AGN can in principle be added as well. 

There are several noteworthy results we obtained with 
our cosmic ray treatment in this study. First of all, our 
simulations of galaxy formation with cosmic ray produc- 
tion by supernovae indicate that cosmic ray pressure can 
play an important role in regulating star formation in 
small galaxies. Here we find a significant reduction of 
the star formation rate compared to the one without CR 
physics, provided cosmic ray production efficiencies of sev- 
eral tens of percent are assumed. In small galaxies, the 
mean densities reached in the ISM stay sufficiently low 
such that the CR pressure can exceed the effective pres- 
sure produced by the thermal supernova feedback. Once 
this occurs, the gas of the ISM is puffed up, quenching the 
star formation rate. Due to the comparatively long cosmic 
ray dissipation timescale, the CR-pressure survives for a 
sufficiently long time in these systems and develops a siz- 



able impact on the star formation rate. In massive galax- 
ies on the other hand, the ISM becomes so dense that the 
CR-pressure is unable to exceed the effective pressure pre- 
dicted by the multi-phase model of ISprineel fc Hernauistl 
ll2003al) .' such that the star formation rates are not altered. 

This effect on star formation also manifests itself in a 
reduction of the cosmic star formation rate density in cos- 
mological simulations of galaxy formation. Here the SFR 
history is reduced at high redshift, where the bulk of star 
formation is dominated by small dwarf galaxies. As the 
star formation shifts to the scale of more massive halos 
towards low redshift, the reduction becomes progressively 
smaller. An interesting implication of the strong effect of 
CR feedback on small galaxies is that this reduces the 
faint-end slope of the resulting galaxy luminosity function, 
an area that continues to be a problematic issue for hy- 
drodynamical simulations of galaxy formation within the 
ACDM scenario. We have indeed detected this flattening, 
although with a weak strength overall. Another tantaliz- 
ing effect of cosmic rays is that they help to keep gas in 
small galaxies more diffuse. This should in principle help 
to alleviate the "angular momentum problem" , which de- 
scribes the problem of an efficient angular momentum loss 
of gas to the dark matter caused by the early collapse of 
large amounts of gas in small halos. It is believed that this 
is a primary reason why present simulations generally fail 
to produce large spiral galaxies at low redshift. Cosmic 
rays physics might help to resolve this problem. 

In simulations where we included cosmic ray injection 
by structure formation shocks, we find that they are pro- 
duced efficiently at high redshift when structure formation 
ensues, driven by the high Mach-number shocks found at 
low to moderate overdensities. At low redshift on the other 
hand, most of the energy is thermalized in weaker shocks 
where the injection efficiency is much smaller. As a re- 
sult, the mean energy content in cosmic rays can reach 
above 40% at redshifts z ~ 5, but drops to ~ 10% at low 
redshift. However, the relative energy content also shows a 
strong density dependence. It is highest at low to moderate 
overdensities and declines continously with density, such 
that deep inside halos, only comparatively little cosmic 
ray energy produced by shock waves survive. An impor- 
tant factor in this trend is the strong density dependence 
of cosmic ray loss processes, and the softer adiabatic index 
of CRs. 

When full cosmological simulations of the formation of 
galaxy clusters are considered, it is therefore not very sur- 
prising that we find that structure formation shocks can 
build up only a comparatively small cosmic ray pressure 
contribution inside clusters. Even at the virial radius this 
contribution reaches only about 10%, but lies much lower 
in the inner parts of the cluster. When radiative cooling 
and cosmic ray production by supernovae are included, we 
find that supernovae can boost the mean CR energy den- 
sity in the cluster, but the averaged contribution still stays 
at the percent level throughout the cluster volume, except 
at places where the gas rapidly cools. Here the CR pres- 
sure can temporarily dominate the pressure and delay the 
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collapse briefly. Nevertheless, we find that CR production 
by supcrnovae and structure formation shocks is unable 
to reduce central cluster cooling flows. Instead, we in fact 
detect a slight increase of the cooling in the 10 14 /i _1 M Q 
cluster we have simulated. This can be understood as a 
result of the higher compressibility of the cluster gas in 
the cosmic ray simulations, leading to an increased cen- 
tral concentration of the gas and an elevated baryon frac- 
tion in the cluster, and thereby to higher cooling of gas 
in the centre overall. Note however that the currently dis- 
cussed AGN feedback for cooling-flow quenching was not 
included in our simulations. On the other hand, the bulk 
of the cluster galaxies experience a reduction of their star 
formation rate when CR feedback is included, such that 
the cluster galaxy luminosity function is expected to de- 
velop a shallower faint-end slope. 

Overall, our results suggest that cosmic ray physics 
is unlikely to drastically modify the physics of galaxy for- 
mation in the ACDM model. However, cosmic rays help in 
areas where current model-building faces important prob- 
lems, like for the faint-end slope of the galaxy luminos- 
ity function and the angular momentum problem. Our 
formalism for treating CRs in cosmological simulations 
should therefore be very valuable for future studies on the 
role of cosmic rays in cosmological structure formation. In 
particular, it would be highly interesting to examine the 
effects of CRs on the metal distribution of the universe, 
or on the dynamics of buoyant bubbles inflated by AGN 
in clusters of galaxies. It will also be important to provide 
an in-depth analysis of the role of cosmic ray diffusion in 
future work. 
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